﻿//201911 大改 符合多种图像格式
//20201016 修改bmp文件读写函数，使不依赖<windows.h>头文件
//20201023 修改图像数据为向量方式存储,使用模板形式重写PImage

#pragma once

//#include "PBaseCommon.h"
#include "PRegion.h"
#include "PString.h"
#include <tchar.h>
#include <typeinfo>

void PImageInitResource();
void PImageExitResource();
void PImageBuildRotateTable();
void PImageDeleteRotateTable();
const byte* PImageRotatePara(int win, int hin, double angle, int& wout, int& hout, PRegion& reg);

inline void type_copy(const byte t1, byte& t2)
{
	t2 = t1;
}
inline void type_copy(const byte t1, char& t2)
{
	t2 = t1;
}
inline void type_copy(const byte t1,uint16_t& t2)
{
	t2 = t1 * 256;
}
inline void type_copy(const byte t1, int16_t& t2)
{
	t2 = t1;
}
inline void type_copy(const byte t1, uint32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const byte t1, int32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const byte t1, float& t2)
{
	t2 = float(t1);
}
inline void type_copy(const byte t1, double& t2)
{
	t2 = double(t1);
}
inline void type_copy(const char t1, byte& t2)
{
	t2 = t1;
}
inline void type_copy(const char t1, char& t2)
{
	t2 = t1;
}
inline void type_copy(const char t1, uint16_t& t2)
{
	t2 = t1;
}
inline void type_copy(const char t1, int16_t& t2)
{
	t2 = t1 * 256;
}
inline void type_copy(const char t1, uint32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const char t1, int32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const char t1, float& t2)
{
	t2 = float(t1);
}
inline void type_copy(const char t1, double& t2)
{
	t2 = double(t1);
}
inline void type_copy(const uint16_t t1, byte& t2)
{
	t2 = byte(t1/256);
}
inline void type_copy(const uint16_t t1, char& t2)
{
	t2 = t1 / 256;
}
inline void type_copy(const uint16_t t1, uint16_t& t2)
{
	t2 = t1;
}
inline void type_copy(const uint16_t t1, int16_t& t2)
{
	t2 = t1;
}
inline void type_copy(const uint16_t t1, uint32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const uint16_t t1, int32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const uint16_t t1, float& t2)
{
	t2 = float(t1);
}
inline void type_copy(const uint16_t t1, double& t2)
{
	t2 = double(t1);
}
inline void type_copy(const int16_t t1, byte& t2)
{
	t2 = byte(uint16_t(t1) / 256);
}
inline void type_copy(const int16_t t1, char& t2)
{
	t2 = char(t1 / 256);
}
inline void type_copy(const int16_t t1, uint16_t& t2)
{
	t2 = t1;
}
inline void type_copy(const int16_t t1, int16_t& t2)
{
	t2 = t1;
}
inline void type_copy(const int16_t t1, uint32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const int16_t t1, int32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const int16_t t1, float& t2)
{
	t2 = float(t1);
}
inline void type_copy(const int16_t t1, double& t2)
{
	t2 = double(t1);
}
inline void type_copy(const uint32_t t1, byte& t2)
{
	if (t1 > 255) t2 = 255;
	else t2 = byte(t1);
}
inline void type_copy(const uint32_t t1, char& t2)
{
	if (t1 > 127) t2 = 127;
	else t2 = char(t1);
}
inline void type_copy(const uint32_t t1, uint16_t& t2)
{
	if (t1 > 65535) t2 = 65535;
	else t2 = uint16_t(t1);
}
inline void type_copy(const uint32_t t1, int16_t& t2)
{
	if (t1 > 32767) t2 = 32767;
	else t2 = int16_t(t1);
}
inline void type_copy(const uint32_t t1, uint32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const uint32_t t1, int32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const uint32_t t1, float& t2)
{
	t2 = float(t1);
}
inline void type_copy(const uint32_t t1, double& t2)
{
	t2 = double(t1);
}
inline void type_copy(const int32_t t1, byte& t2)
{
	if (t1 > 255) t2 = 255;
	else if (t1 < 0) t2 = 0;
	else t2 = byte(t1);
}
inline void type_copy(const int32_t t1, char& t2)
{
	if (t1 > 127) t2 = 127;
	else if (t1 < -128) t2 = -128;
	else t2 = char(t1);
}
inline void type_copy(const int32_t t1, uint16_t& t2)
{
	if (t1 > 65535) t2 = 65535;
	else if (t1 < 0) t2 = 0;
	else t2 = uint16_t(t1);
}
inline void type_copy(const int32_t t1, int16_t& t2)
{
	if (t1 > 32767) t2 = 32767;
	else if (t1 < -32768) t2 = -32768;
	else t2 = int16_t(t1);
}
inline void type_copy(const int32_t t1, uint32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const int32_t t1, int32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const int32_t t1, float& t2)
{
	t2 = float(t1);
}
inline void type_copy(const int32_t t1, double& t2)
{
	t2 = double(t1);
}
inline void type_copy(const float t1, byte& t2)
{
	if (t1 > 255.0) t2 = 255;
	else if (t1 < 0.0) t2 = 0;
	else t2 = byte(t1);
}
inline void type_copy(const float t1, char& t2)
{
	if (t1 > 127.0) t2 = 127;
	else if (t1 < -128.0) t2 = -128;
	else t2 = char(t1);
}
inline void type_copy(const float t1, uint16_t& t2)
{
	if (t1 > 65535.0) t2 = 65535;
	else if (t1 < 0.0) t2 = 0;
	else t2 = uint16_t(t1);
}
inline void type_copy(const float t1, int16_t& t2)
{
	if (t1 > 32767.0) t2 = 32767;
	else if (t1 < -32768.0) t2 = -32768;
	else t2 = int16_t(t1);
}
inline void type_copy(const float t1, uint32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const float t1, int32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const float t1, float& t2)
{
	t2 = t1;
}
inline void type_copy(const float t1, double& t2)
{
	t2 = t1;
}
inline void type_copy(const double t1, byte& t2)
{
	if (t1 > 255.0) t2 = 255;
	else if (t1 < 0.0) t2 = 0;
	else t2 = byte(t1);
}
inline void type_copy(const double t1, char& t2)
{
	if (t1 > 127.0) t2 = 127;
	else if (t1 < -128.0) t2 = -128;
	else t2 = char(t1);
}
inline void type_copy(const double t1, uint16_t& t2)
{
	if (t1 > 65535.0) t2 = 65535;
	else if (t1 < 0.0) t2 = 0;
	else t2 = uint16_t(t1);
}
inline void type_copy(const double t1, int16_t& t2)
{
	if (t1 > 32767.0) t2 = 32767;
	else if (t1 < -32768.0) t2 = -32768;
	else t2 = int16_t(t1);
}
inline void type_copy(const double t1, uint32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const double t1, int32_t& t2)
{
	t2 = t1;
}
inline void type_copy(const double t1, float& t2)
{
	t2 = t1;
}
inline void type_copy(const double t1, double& t2)
{
	t2 = t1;
}

template<typename T> class ptv1
{
public:
	typedef PDSum1 DSUM;
	typedef T ITEM;
	ptv1()
	{
		v = 0;
	}
	ptv1(T v1)
	{
		v = v1;
	}
	ptv1(const pstring& value)
	{
		v=value.toInt();
	}
	ptv1(const ptv1& t)
	{
		v = t.v;
	}
	static inline int channels()
	{
		return 1;
	}
	static inline int size()
	{
		return sizeof(T);
	}
	static inline const type_info& type()
	{
		return typeid(T);
	}
	inline T& operator[](int i)
	{
		return v;
	}
	inline T operator[](int i) const
	{
		return v;
	}
	inline void zero()
	{
		v = 0;
	}
	inline void set(T v1)
	{
		v = v1;
	}
	inline void set(const ptv1& t)
	{
		v = t.v;
	}
	inline void put(DSUM& s) const
	{
		s[0] = v;
	}
	inline void put(T* p) const
	{
		*p = v;
	}
	inline operator DSUM() const
	{
		return DSUM(double(v));
	}
	inline DSUM operator+(const ptv1& t) const
	{
		return DSUM(double(v) + t.v);
	}
	inline DSUM operator-(const ptv1& t) const
	{
		return DSUM(double(v) - t.v);
	}
	inline DSUM operator*(const ptv1& t) const
	{
		return DSUM(double(v) * t.v);
	}
	inline DSUM operator/(const ptv1& t) const
	{
		return DSUM(double(v) / t.v);
	}
	inline DSUM operator+(double d) const
	{
		return DSUM(double(v) + d);
	}
	inline DSUM operator-(double d) const
	{
		return DSUM(double(v) - d);
	}
	inline DSUM operator*(double d) const
	{
		return DSUM(double(v) * d);
	}
	inline DSUM operator/(double d) const
	{
		return DSUM(double(v) / d);
	}
	inline ptv1& operator+=(const ptv1& t)
	{
		type_copy(v + t.v, v);
		return *this;
	}
	inline ptv1& operator-=(const ptv1& t)
	{
		type_copy(v - t.v, v);
		return *this;
	}
	inline ptv1& operator*=(const ptv1& t)
	{
		type_copy(v * t.v, v);
		return *this;
	}
	inline ptv1& operator/=(const ptv1& t)
	{
		type_copy(v / t.v, v);
		return *this;
	}
	inline ptv1& operator+=(double d)
	{
		type_copy(v + d, v);
		return *this;
	}
	inline ptv1& operator-=(double d)
	{
		type_copy(v - d, v);
		return *this;
	}
	inline ptv1& operator*=(double d)
	{
		type_copy(v * d, v);
		return *this;
	}
	inline ptv1& operator/=(double d)
	{
		type_copy(v / d, v);
		return *this;
	}
	inline ptv1 operator&(const ptv1& t) const
	{
		return ptv1(v & t.v);
	}
	inline ptv1 operator|(const ptv1& t) const
	{
		return ptv1(v | t.v);
	}
	inline ptv1 operator~() const
	{
		return ptv1(~v);
	}
	inline ptv1& operator&=(const ptv1& t)
	{
		v &= t.v;
		return *this;
	}
	inline ptv1& operator|=(const ptv1& t)
	{
		v |= t.v;
		return *this;
	}
	inline ptv1& operator=(const DSUM& s)
	{
		type_copy(s[0], v);
		return *this;
	}	
	inline ptv1& operator=(const ptv1& t)
	{
		v = t.v;
		return *this;
	}
	inline bool operator==(const ptv1& t) const
	{
		return (v==t.v);
	}
	inline bool operator<(const ptv1& t) const
	{
		return (v<t.v);
	}
	inline bool operator<=(const ptv1& t) const
	{
		return (v<=t.v);
	}
	inline bool operator>(const ptv1& t) const
	{
		return (v>t.v);
	}
	inline bool operator>=(const ptv1& t) const
	{
		return (v>=t.v);
	}
	inline bool operator==(const DSUM& t) const
	{
		return (t.v==v);
	}
	inline bool operator<(const DSUM& t) const
	{
		return (t.v>v);
	}
	inline bool operator<=(const DSUM& t) const
	{
		return (t.v>=v);
	}
	inline bool operator>(const DSUM& t) const
	{
		return (t.v<v);
	}
	inline bool operator>=(const DSUM& t) const
	{
		return (t.v<=v);
	}
	inline double distance(const ptv1 &t) const
	{
		double s(double(v) - t.v);
		return sqrt(s*s);
	}
	inline double distance2(const ptv1 &t) const
	{
		double s(double(v) - t.v);
		return s*s;
	}
	inline double average() const
	{
		return double(v);
	}
	inline T vmin() const
	{
		return v;
	}
	inline T vmax() const
	{
		return v;
	}
	inline DSUM gamma(double c,double g) const
	{
		return DSUM(c*pow(double(v),g));
	}
	inline DSUM clog(double c) const
	{
		return DSUM(c * log(1.0+v));
	}
	T v;
};

template<typename T> class ptv2
{
public:
	typedef PDSum2 DSUM;
	typedef T ITEM;
	ptv2()
	{
		v[0] = 0;
		v[1] = 0;
	}
	ptv2(T v1,T v2)
	{
		v[0] = v1;
		v[1] = v2;
	}
	ptv2(const ptv2& t)
	{
		v[0] = t[0];
		v[1] = t[1];
	}
	ptv2(const pstring& value)
	{
		int i = 0;
		v[1] = value.getRange(i, ',').toInt();
		v[0] = value.getRange(i, ',').toInt();
	}
	static inline int channels()
	{
		return 2;
	}
	static inline int size()
	{
		return sizeof(T);
	}
	static inline const type_info& type()
	{
		return typeid(T);
	}
	inline T& operator[](int i)
	{
		return v[i];
	}
	inline T operator[](int i) const
	{
		return v[i];
	}
	inline void zero()
	{
		v[0] = 0;
		v[1] = 0;
	}
	inline void set(T v1,T v2)
	{
		v[0] = v1;
		v[1] = v2;
	}
	inline void set(const ptv2& t)
	{
		v[0] = t[0];
		v[1] = t[1];
	}
	inline void put(DSUM& s) const
	{
		s[0] = v[0];
		s[1] = v[1];
	}
	inline void put(T* p) const
	{
		p[0] = v[0];
		p[1] = v[1];
	}
	inline operator DSUM() const
	{
		return DSUM(double(v[0]),double(v[1]));
	}
	inline DSUM operator+(const ptv2& t) const
	{
		return DSUM(double(v[0]) + t[0],double(v[1]) + t[1]);
	}
	inline DSUM operator-(const ptv2& t) const
	{
		return DSUM(double(v[0]) - t[0],double(v[1]) - t[1]);
	}
	inline DSUM operator*(const ptv2& t) const
	{
		return DSUM(double(v[0]) * t[0],double(v[1]) * t[1]);
	}
	inline DSUM operator/(const ptv2& t) const
	{
		return DSUM(double(v[0]) / t[0],double(v[1]) / t[1]);
	}
	inline DSUM operator+(double d) const
	{
		return DSUM(double(v[0]) + d,double(v[1]) + d);
	}
	inline DSUM operator-(double d) const
	{
		return DSUM(double(v[0]) - d,double(v[1]) - d);
	}
	inline DSUM operator*(double d) const
	{
		return DSUM(double(v[0]) * d,double(v[1]) * d);
	}
	inline DSUM operator/(double d) const
	{
		return DSUM(double(v[0]) / d,double(v[1]) / d);
	}
	inline ptv2& operator+=(const ptv2& t)
	{
		type_copy(v[0] + t[0], v[0]);
		type_copy(v[1] + t[1], v[1]);
		return *this;
	}
	inline ptv2& operator-=(const ptv2& t)
	{
		type_copy(v[0] - t[0], v[0]);
		type_copy(v[1] - t[1], v[1]);
		return *this;
	}
	inline ptv2& operator*=(const ptv2& t)
	{
		type_copy(v[0] * t[0], v[0]);
		type_copy(v[1] * t[1], v[1]);
		return *this;
	}
	inline ptv2& operator/=(const ptv2& t)
	{
		type_copy(v[0] / t[0], v[0]);
		type_copy(v[1] / t[1], v[1]);
		return *this;
	}
	inline ptv2& operator+=(double d)
	{
		type_copy(v[0] + d, v[0]);
		type_copy(v[1] + d, v[1]);
		return *this;
	}
	inline ptv2& operator-=(double d)
	{
		type_copy(v[0] - d, v[0]);
		type_copy(v[1] - d, v[1]);
		return *this;
	}
	inline ptv2& operator*=(double d)
	{
		type_copy(v[0] * d, v[0]);
		type_copy(v[1] * d, v[1]);
		return *this;
	}
	inline ptv2& operator/=(double d)
	{
		type_copy(v[0] / d, v[0]);
		type_copy(v[1] / d, v[1]);
		return *this;
	}
	inline ptv2 operator&(const ptv2& t) const
	{
		return ptv2(v[0] & t[0],v[1] & t[1]);
	}
	inline ptv2 operator|(const ptv2& t) const
	{
		return ptv2(v[0] | t[0],v[1] | t[1]);
	}
	inline ptv2 operator~() const
	{
		return ptv2(~v[0],~v[1]);
	}
	inline ptv2& operator&=(const ptv2& t)
	{
		v[0] &= t[0];
		v[1] &= t[1];
		return *this;
	}
	inline ptv2& operator|=(const ptv2& t)
	{
		v[0] |= t[0];
		v[1] |= t[1];
		return *this;
	}
	inline ptv2& operator=(const DSUM& s)
	{
		type_copy(s[0], v[0]);
		type_copy(s[1], v[1]);
		return *this;
	}	
	inline ptv2& operator=(const ptv2& t)
	{
		v[0] = t[0];
		v[1] = t[1];
		return *this;
	}
	inline bool operator==(const ptv2& t) const
	{
		return (v[0]==t[0] && v[1]==t[1]);
	}
	inline bool operator<(const ptv2& t) const
	{
		return (v[0]<t[0] && v[1]<t[1]);
	}
	inline bool operator<=(const ptv2& t) const
	{
		return (v[0]<=t[0] && v[1]<=t[1]);
	}
	inline bool operator>(const ptv2& t) const
	{
		return (v[0]>t[0] && v[1]>t[1]);
	}
	inline bool operator>=(const ptv2& t) const
	{
		return (v[0]>=t[0] && v[1]>=t[1]);
	}
	inline bool operator==(const DSUM& t) const
	{
		return (t[0]==v[0] && t[1]==v[1]);
	}
	inline bool operator<(const DSUM& t) const
	{
		return (t[0]>v[0] && t[1]>v[1]);
	}
	inline bool operator<=(const DSUM& t) const
	{
		return (t[0]>=v[0] && t[1]>=v[1]);
	}
	inline bool operator>(const DSUM& t) const
	{
		return (t[0]<v[0] && t[1]<v[1]);
	}
	inline bool operator>=(const DSUM& t) const
	{
		return (t[0]<=v[0] && t[1]<=v[1]);
	}
	inline double distance(const ptv2 &t) const
	{
		DSUM s(double(v[0]) - t[0],double(v[1]) - t[1]);
		return sqrt(s[0]*s[0]+s[1]*s[1]);
	}
	inline double distance2(const ptv2 &t) const
	{
		DSUM s(double(v[0]) - t[0],double(v[1]) - t[1]);
		return s[0]*s[0]+s[1]*s[1];
	}
	inline double average() const
	{
		return (double(v[0])+v[1])/2.0;
	}
	inline T vmin() const
	{
		if (v[0] <= v[1]) return v[0];
		else return v[1];
	}
	inline T vmax() const
	{
		if (v[0] >= v[1]) return v[0];
		else return v[1];
	}
	inline DSUM gamma(double c,double g) const
	{
		return DSUM(c*pow(double(v[0]),g),c*pow(double(v[1]),g));
	}
	inline DSUM clog(double c) const
	{
		return DSUM(c * log(1.0 + v[0]), c * log(1.0 + v[1]));
	}
	T v[2];
};

template<typename T> class ptv3
{
public:
	typedef PDSum3 DSUM;
	typedef T ITEM;
	ptv3()
	{
		v[0] = 0;
		v[1] = 0;
		v[2] = 0;
	}
	ptv3(T v1,T v2,T v3)
	{
		v[0] = v1;
		v[1] = v2;
		v[2] = v3;
	}
	ptv3(const ptv3& t)
	{
		v[0] = t[0];
		v[1] = t[1];
		v[2] = t[2];
	}
	ptv3(const pstring& value)
	{
		int i = 0;
		v[2] = value.getRange(i, ',').toInt();
		v[1] = value.getRange(i, ',').toInt();
		v[0] = value.getRange(i, ',').toInt();
	}
	static inline int channels()
	{
		return 3;
	}
	static inline int size()
	{
		return sizeof(T);
	}
	static inline const type_info& type()
	{
		return typeid(T);
	}
	inline T& operator[](int i)
	{
		return v[i];
	}
	inline T operator[](int i) const
	{
		return v[i];
	}
	inline void zero()
	{
		v[0] = 0;
		v[1] = 0;
		v[2] = 0;
	}
	inline void set(T v1,T v2,T v3)
	{
		v[0] = v1;
		v[1] = v2;
		v[2] = v3;
	}
	inline void set(const ptv3& t)
	{
		v[0] = t[0];
		v[1] = t[1];
		v[2] = t[2];
	}
	inline void put(DSUM& s) const
	{
		s[0] = v[0];
		s[1] = v[1];
		s[2] = v[2];
	}
	inline void put(T* p) const
	{
		p[0] = v[0];
		p[1] = v[1];
		p[2] = v[2];
	}
	inline operator DSUM() const
	{
		return DSUM(double(v[0]),double(v[1]),double(v[2]));
	}
	inline DSUM operator+(const ptv3& t) const
	{
		return DSUM(double(v[0]) + t[0],double(v[1]) + t[1],double(v[2]) + t[2]);
	}
	inline DSUM operator-(const ptv3& t) const
	{
		return DSUM(double(v[0]) - t[0],double(v[1]) - t[1],double(v[2]) - t[2]);
	}
	inline DSUM operator*(const ptv3& t) const
	{
		return DSUM(double(v[0]) * t[0],double(v[1]) * t[1],double(v[2]) * t[2]);
	}
	inline DSUM operator/(const ptv3& t) const
	{
		return DSUM(double(v[0]) / t[0],double(v[1]) / t[1],double(v[2]) / t[2]);
	}
	inline DSUM operator+(double d) const
	{
		return DSUM(double(v[0]) + d,double(v[1]) + d,double(v[2]) + d);
	}
	inline DSUM operator-(double d) const
	{
		return DSUM(double(v[0]) - d,double(v[1]) - d,double(v[2]) - d);
	}
	inline DSUM operator*(double d) const
	{
		return DSUM(double(v[0]) * d,double(v[1]) * d,double(v[2]) * d);
	}
	inline DSUM operator/(double d) const
	{
		return DSUM(double(v[0]) / d,double(v[1]) / d,double(v[2]) / d);
	}
	inline ptv3& operator+=(const ptv3& t)
	{
		type_copy(v[0] + t[0], v[0]);
		type_copy(v[1] + t[1], v[1]);
		type_copy(v[2] + t[2], v[2]);
		return *this;
	}
	inline ptv3& operator-=(const ptv3& t)
	{
		type_copy(v[0] - t[0], v[0]);
		type_copy(v[1] - t[1], v[1]);
		type_copy(v[2] - t[2], v[2]);
		return *this;
	}
	inline ptv3& operator*=(const ptv3& t)
	{
		type_copy(v[0] * t[0], v[0]);
		type_copy(v[1] * t[1], v[1]);
		type_copy(v[2] * t[2], v[2]);
		return *this;
	}
	inline ptv3& operator/=(const ptv3& t)
	{
		type_copy(v[0] / t[0], v[0]);
		type_copy(v[1] / t[1], v[1]);
		type_copy(v[2] / t[2], v[2]);
		return *this;
	}
	inline ptv3& operator+=(double d)
	{
		type_copy(v[0] + d, v[0]);
		type_copy(v[1] + d, v[1]);
		type_copy(v[2] + d, v[2]);
		return *this;
	}
	inline ptv3& operator-=(double d)
	{
		type_copy(v[0] - d, v[0]);
		type_copy(v[1] - d, v[1]);
		type_copy(v[2] - d, v[2]);
		return *this;
	}
	inline ptv3& operator*=(double d)
	{
		type_copy(v[0] * d, v[0]);
		type_copy(v[1] * d, v[1]);
		type_copy(v[2] * d, v[2]);
		return *this;
	}
	inline ptv3& operator/=(double d)
	{
		type_copy(v[0] / d, v[0]);
		type_copy(v[1] / d, v[1]);
		type_copy(v[2] / d, v[2]);
		return *this;
	}
	inline ptv3 operator&(const ptv3& t) const
	{
		return ptv3(v[0] & t[0],v[1] & t[1],v[2] & t[2]);
	}
	inline ptv3 operator|(const ptv3& t) const
	{
		return ptv3(v[0] | t[0],v[1] | t[1],v[2] | t[2]);
	}
	inline ptv3 operator~() const
	{
		return ptv3(~v[0],~v[1],~v[2]);
	}
	inline ptv3& operator&=(const ptv3& t)
	{
		v[0] &= t[0];
		v[1] &= t[1];
		v[2] &= t[2];
		return *this;
	}
	inline ptv3& operator|=(const ptv3& t)
	{
		v[0] |= t[0];
		v[1] |= t[1];
		v[2] |= t[2];
		return *this;
	}
	inline ptv3& operator=(const DSUM& s)
	{
		type_copy(s[0], v[0]);
		type_copy(s[1], v[1]);
		type_copy(s[2], v[2]);
		return *this;
	}	
	inline ptv3& operator=(const ptv3& t)
	{
		v[0] = t[0];
		v[1] = t[1];
		v[2] = t[2];
		return *this;
	}
	inline bool operator==(const ptv3& t) const
	{
		return (v[0]==t[0] && v[1]==t[1] && v[2]==t[2]);
	}
	inline bool operator<(const ptv3& t) const
	{
		return (v[0]<t[0] && v[1]<t[1] && v[2]<t[2]);
	}
	inline bool operator<=(const ptv3& t) const
	{
		return (v[0]<=t[0] && v[1]<=t[1] && v[2]<=t[2]);
	}
	inline bool operator>(const ptv3& t) const
	{
		return (v[0]>t[0] && v[1]>t[1] && v[2]>t[2]);
	}
	inline bool operator>=(const ptv3& t) const
	{
		return (v[0]>=t[0] && v[1]>=t[1] && v[2]>=t[2]);
	}
	inline bool operator==(const DSUM& t) const
	{
		return (t[0]==v[0] && t[1]==v[1] && t[2]==v[2]);
	}
	inline bool operator<(const DSUM& t) const
	{
		return (t[0]>v[0] && t[1]>v[1] && t[2]>v[2]);
	}
	inline bool operator<=(const DSUM& t) const
	{
		return (t[0]>=v[0] && t[1]>=v[1] && t[2]>=v[2]);
	}
	inline bool operator>(const DSUM& t) const
	{
		return (t[0]<v[0] && t[1]<v[1] && t[2]<v[2]);
	}
	inline bool operator>=(const DSUM& t) const
	{
		return (t[0]<=v[0] && t[1]<=v[1] && t[2]<=v[2]);
	}
	inline double distance(const ptv3 &t) const
	{
		DSUM s(double(v[0]) - t[0],double(v[1]) - t[1],double(v[2]) - t[2]);
		return sqrt(s[0]*s[0]+s[1]*s[1]+s[2]*s[2]);
	}
	inline double distance2(const ptv3 &t) const
	{
		DSUM s(double(v[0]) - t[0],double(v[1]) - t[1],double(v[2]) - t[2]);
		return s[0]*s[0]+s[1]*s[1]+s[2]*s[2];
	}
	inline double average() const
	{
		return (double(v[0])+v[1]+v[2])/3.0;
	}
	inline T vmin() const
	{
		if (v[0] <= v[1] && v[0] <= v[2]) return v[0];
		else if (v[1] <= v[2]) return v[1];
		else return v[2];
	}
	inline T vmax() const
	{
		if (v[0] >= v[1] && v[0] >= v[2]) return v[0];
		else if (v[1] >= v[2]) return v[1];
		else return v[2];
	}
	inline DSUM gamma(double c,double g) const
	{
		return DSUM(c*pow(double(v[0]),g),c*pow(double(v[1]),g),c*pow(double(v[2]),g));
	}
	inline DSUM clog(double c) const
	{
		return DSUM(c * log(1.0 + v[0]), c * log(1.0 + v[1]), c * log(1.0 + v[2]));
	}
	inline void hsi(byte& h, byte& s, byte& i) const
	{
		double r = v[2], g = v[1], b = v[0];
		double dh = acos(0.5 * ((r - g) + (r - b)) / sqrt((r - g) * (r - g) + (r - b) * (g - b)));
		if (b > g) dh = pi * 2.0 - dh;
		double m = min(r, g);
		m = min(m, b);
		double ds = 1.0 - 3.0 * m / (r + g + b);
		double di = (r + g + b) / 3.0;
		s = byte(ds * 255.5);
		if (s == 0) h = 0;
		else h = byte(dh / (pi * 2) * 255.5);
		i= byte(di);
	}
	inline void rgb(byte& r, byte& g, byte& b)
	{
		double h = v[0] * pi * 2.0 / 255.0, s = v[1]/255.0, i = v[2];
		double dr, dg, db;
	
		if (h < 2.0*pi/3.0)
		{
			db = i * (1.0 - s);
			dr = i * (1.0 + s * cos(h) / cos(pi / 3.0 - h));
			dg = i * 3.0 - (db + dr);
		}
		else if (h < 4.0 * pi / 3.0)
		{
			h -= 2.0 * pi / 3.0;
			dr = i * (1.0 - s);
			dg = i * (1.0 + s * cos(h) / cos(pi / 3.0 - h));
			db = i * 3.0 - (dg + dr);
		}
		else
		{
			h -= 4.0 * pi / 3.0;
			dg = i * (1.0 - s);
			db = i * (1.0 + s * cos(h) / cos(pi / 3.0 - h));
			dr = i * 3.0 - (dg + db);
		}
		type_copy(dr, r);
		type_copy(dg, g);
		type_copy(db, b);
	}
	T v[3];
};

template<typename T> class ptv4
{
public:
	typedef PDSum4 DSUM;
	typedef T ITEM;
	ptv4()
	{
		v[0] = 0;
		v[1] = 0;
		v[2] = 0;
		v[3] = 0;
	}
	ptv4(T v1,T v2,T v3,T v4)
	{
		v[0] = v1;
		v[1] = v2;
		v[2] = v3;
		v[3] = v4;
	}
	ptv4(const ptv4& t)
	{
		v[0] = t[0];
		v[1] = t[1];
		v[2] = t[2];
		v[3] = t[3];
	}
	ptv4(const pstring& value)
	{
		int i = 0;
		v[3] = value.getRange(i, ',').toInt();
		v[2] = value.getRange(i, ',').toInt();
		v[1] = value.getRange(i, ',').toInt();
		v[0] = value.getRange(i, ',').toInt();
	}
	static inline int channels()
	{
		return 4;
	}
	static inline int size()
	{
		return sizeof(T);
	}
	static inline const type_info& type()
	{
		return typeid(T);
	}
	inline T& operator[](int i)
	{
		return v[i];
	}
	inline T operator[](int i) const
	{
		return v[i];
	}
	inline void zero()
	{
		v[0] = 0;
		v[1] = 0;
		v[2] = 0;
		v[3] = 0;
	}
	inline void set(T v1,T v2,T v3,T v4)
	{
		v[0] = v1;
		v[1] = v2;
		v[2] = v3;
		v[3] = v4;
	}
	inline void set(const ptv4& t)
	{
		v[0] = t[0];
		v[1] = t[1];
		v[2] = t[2];
		v[3] = t[3];
	}
	inline void put(DSUM& s) const
	{
		s[0] = v[0];
		s[1] = v[1];
		s[2] = v[2];
		s[3] = v[3];
	}
	inline void put(T* p) const
	{
		p[0] = v[0];
		p[1] = v[1];
		p[2] = v[2];
		p[3] = v[3];
	}
	inline operator DSUM() const
	{
		return DSUM(double(v[0]),double(v[1]),double(v[2]),double(v[3]));
	}
	inline DSUM operator+(const ptv4& t) const
	{
		return DSUM(double(v[0]) + t[0],double(v[1]) + t[1],double(v[2]) + t[2],double(v[3]) + t[3]);
	}
	inline DSUM operator-(const ptv4& t) const
	{
		return DSUM(double(v[0]) - t[0],double(v[1]) - t[1],double(v[2]) - t[2],double(v[3]) - t[3]);
	}
	inline DSUM operator*(const ptv4& t) const
	{
		return DSUM(double(v[0]) * t[0],double(v[1]) * t[1],double(v[2]) * t[2],double(v[3]) * t[3]);
	}
	inline DSUM operator/(const ptv4& t) const
	{
		return DSUM(double(v[0]) / t[0],double(v[1]) / t[1],double(v[2]) / t[2],double(v[3]) / t[3]);
	}
	inline DSUM operator+(double d) const
	{
		return DSUM(double(v[0]) + d,double(v[1]) + d,double(v[2]) + d,double(v[3]) + d);
	}
	inline DSUM operator-(double d) const
	{
		return DSUM(double(v[0]) - d,double(v[1]) - d,double(v[2]) - d,double(v[3]) - d);
	}
	inline DSUM operator*(double d) const
	{
		return DSUM(double(v[0]) * d,double(v[1]) * d,double(v[2]) * d,double(v[3]) * d);
	}
	inline DSUM operator/(double d) const
	{
		return DSUM(double(v[0]) / d,double(v[1]) / d,double(v[2]) / d,double(v[3]) / d);
	}
	inline ptv4& operator+=(const ptv4& t)
	{
		type_copy(v[0] + t[0], v[0]);
		type_copy(v[1] + t[1], v[1]);
		type_copy(v[2] + t[2], v[2]);
		type_copy(v[3] + t[3], v[3]);
		return *this;
	}
	inline ptv4& operator-=(const ptv4& t)
	{
		type_copy(v[0] - t[0], v[0]);
		type_copy(v[1] - t[1], v[1]);
		type_copy(v[2] - t[2], v[2]);
		type_copy(v[3] - t[3], v[3]);
		return *this;
	}
	inline ptv4& operator*=(const ptv4& t)
	{
		type_copy(v[0] * t[0], v[0]);
		type_copy(v[1] * t[1], v[1]);
		type_copy(v[2] * t[2], v[2]);
		type_copy(v[3] * t[3], v[3]);
		return *this;
	}
	inline ptv4& operator/=(const ptv4& t)
	{
		type_copy(v[0] / t[0], v[0]);
		type_copy(v[1] / t[1], v[1]);
		type_copy(v[2] / t[2], v[2]);
		type_copy(v[3] / t[3], v[3]);
		return *this;
	}
	inline ptv4& operator+=(double d)
	{
		type_copy(v[0] + d, v[0]);
		type_copy(v[1] + d, v[1]);
		type_copy(v[2] + d, v[2]);
		type_copy(v[3] + d, v[3]);
		return *this;
	}
	inline ptv4& operator-=(double d)
	{
		type_copy(v[0] - d, v[0]);
		type_copy(v[1] - d, v[1]);
		type_copy(v[2] - d, v[2]);
		type_copy(v[3] - d, v[3]);
		return *this;
	}
	inline ptv4& operator*=(double d)
	{
		type_copy(v[0] * d, v[0]);
		type_copy(v[1] * d, v[1]);
		type_copy(v[2] * d, v[2]);
		type_copy(v[3] * d, v[3]);
		return *this;
	}
	inline ptv4& operator/=(double d)
	{
		type_copy(v[0] / d, v[0]);
		type_copy(v[1] / d, v[1]);
		type_copy(v[2] / d, v[2]);
		type_copy(v[3] / d, v[3]);
		return *this;
	}
	inline ptv4 operator&(const ptv4& t) const
	{
		return ptv4(v[0] & t[0],v[1] & t[1],v[2] & t[2],v[3] & t[3]);
	}
	inline ptv4 operator|(const ptv4& t) const
	{
		return ptv4(v[0] | t[0],v[1] | t[1],v[2] | t[2],v[3] | t[3]);
	}
	inline ptv4 operator~() const
	{
		return ptv3(~v[0],~v[1],~v[2],~v[3]);
	}
	inline ptv4& operator&=(const ptv4& t)
	{
		v[0] &= t[0];
		v[1] &= t[1];
		v[2] &= t[2];
		v[3] &= t[3];
		return *this;
	}
	inline ptv4& operator|=(const ptv4& t)
	{
		v[0] |= t[0];
		v[1] |= t[1];
		v[2] |= t[2];
		v[3] |= t[3];
		return *this;
	}
	inline ptv4& operator=(const DSUM& s)
	{
		type_copy(s[0], v[0]);
		type_copy(s[1], v[1]);
		type_copy(s[2], v[2]);
		type_copy(s[3], v[3]);
		return *this;
	}	
	inline ptv4& operator=(const ptv4& t)
	{
		v[0] = t[0];
		v[1] = t[1];
		v[2] = t[2];
		v[3] = t[3];
		return *this;
	}
	inline bool operator==(const ptv4& t) const
	{
		return (v[0]==t[0] && v[1]==t[1] && v[2]==t[2] && v[3]==t[3]);
	}
	inline bool operator<(const ptv4& t) const
	{
		return (v[0]<t[0] && v[1]<t[1] && v[2]<t[2] && v[3]<t[3]);
	}
	inline bool operator<=(const ptv4& t) const
	{
		return (v[0]<=t[0] && v[1]<=t[1] && v[2]<=t[2] && v[3]<=t[3]);
	}
	inline bool operator>(const ptv4& t) const
	{
		return (v[0]>t[0] && v[1]>t[1] && v[2]>t[2] && v[3]>t[3]);
	}
	inline bool operator>=(const ptv4& t) const
	{
		return (v[0]>=t[0] && v[1]>=t[1] && v[2]>=t[2] && v[3]>=t[3]);
	}
	inline bool operator==(const DSUM& t) const
	{
		return (t[0]==v[0] && t[1]==v[1] && t[2]==v[2] && t[3]==v[3]);
	}
	inline bool operator<(const DSUM& t) const
	{
		return (t[0]>v[0] && t[1]>v[1] && t[2]>v[2] && t[3]>v[3]);
	}
	inline bool operator<=(const DSUM& t) const
	{
		return (t[0]>=v[0] && t[1]>=v[1] && t[2]>=v[2] && t[3]>=v[3]);
	}
	inline bool operator>(const DSUM& t) const
	{
		return (t[0]<v[0] && t[1]<v[1] && t[2]<v[2] && t[3]<v[3]);
	}
	inline bool operator>=(const DSUM& t) const
	{
		return (t[0]<=v[0] && t[1]<=v[1] && t[2]<=v[2] && t[3]<=v[3]);
	}
	inline double distance(const ptv4 &t) const
	{
		DSUM s(double(v[0]) - t[0],double(v[1]) - t[1],
		double(v[2]) - t[2],double(v[3]) - t[3]);
		return sqrt(s[0]*s[0]+s[1]*s[1]+s[2]*s[2]+s[3]*s[3]);
	}
	inline double distance2(const ptv4 &t) const
	{
		DSUM s(double(v[0]) - t[0],double(v[1]) - t[1],
		double(v[2]) - t[2],double(v[3]) - t[3]);
		return s[0]*s[0]+s[1]*s[1]+s[2]*s[2]+s[3]*s[3];
	}
	inline double average() const
	{
		return (double(v[0])+v[1]+v[2]+v[3])/4.0;
	}
	inline T vmin() const
	{
		if (v[0] <= v[1] && v[0] <= v[2] && v[0] <= v[3]) return v[0];
		else if (v[1] <= v[2] && v[1] <= v[3]) return v[1];
		else if (v[2] <= v[3]) return v[2];
		else return v[3];
	}
	inline T vmax() const
	{
		if (v[0] >= v[1] && v[0] >= v[2] && v[0] >= v[3]) return v[0];
		else if (v[1] >= v[2] && v[1] >= v[3]) return v[1];
		else if (v[2] >= v[3]) return v[2];
		else return v[3];
	}
	inline DSUM gamma(double c,double g) const
	{
		return DSUM(c*pow(double(v[0]),g),c*pow(double(v[1]),g),c*pow(double(v[2]),g),c*pow(double(v[3]),g));
	}
	inline DSUM clog(double c) const
	{
		return DSUM(c * log(1.0 + v[0]), c * log(1.0 + v[1]), c * log(1.0 + v[2]), c * log(1.0 + v[3]));
	}
	T v[4];
};

template<typename T, int N> class ptvn
{
public:
	typedef PDNSum<N> DSUM;
	typedef T ITEM;
	ptvn()
	{
		for (int i = 0; i < N; ++i) v[i] = 0;
	}
	ptvn(T* vp)
	{
		for (int i = 0; i < N; ++i) v[i] = vp[i];
	}
	ptvn(const ptvn& t)
	{
		for (int i = 0; i < N; ++i) v[i] = t[i];
	}
	ptvn(const pstring& value)
	{
		int index = 0;
		for (int i = N - 1; i >= 0; --i)
		{
			v[i] = value.getRange(index, ',').toInt();
		}
	}
	static inline int channels()
	{
		return N;
	}
	static inline int size()
	{
		return sizeof(T);
	}
	static inline const type_info& type()
	{
		return typeid(T);
	}
	inline T& operator[](int i)
	{
		return v[i];
	}
	inline T operator[](int i) const
	{
		return v[i];
	}
	inline void zero()
	{
		for (int i = 0; i < N; ++i) v[i] = 0;
	}
	inline void set(T* vp)
	{
		for (int i = 0; i < N; ++i) v[i] = vp[i];
	}
	inline void set(const ptvn& t)
	{
		for (int i = 0; i < N; ++i) v[i] = t[i];
	}
	inline void put(DSUM& s) const
	{
		for (int i = 0; i < N; ++i) s[i] = v[i];
	}
	inline void put(T* p) const
	{
		for (int i = 0; i < N; ++i) p[i] = v[i];
	}
	inline operator DSUM() const
	{
		DSUM s;
		for (int i = 0; i < N; ++i) s[i] = v[i];
		return s;
	}
	inline DSUM operator+(const ptvn& t) const
	{
		DSUM s;
		for (int i = 0; i < N; ++i) s[i] = double(v[i]) + t[i];
		return s;
	}
	inline DSUM operator-(const ptvn& t) const
	{
		DSUM s;
		for (int i = 0; i < N; ++i) s[i] = double(v[i]) - t[i];
		return s;
	}
	inline DSUM operator*(const ptvn& t) const
	{
		DSUM s;
		for (int i = 0; i < N; ++i) s[i] = double(v[i]) * t[i];
		return s;
	}
	inline DSUM operator/(const ptvn& t) const
	{
		DSUM s;
		for (int i = 0; i < N; ++i) s[i] = double(v[i]) / t[i];
		return s;
	}
	inline DSUM operator+(double d) const
	{
		DSUM s;
		for (int i = 0; i < N; ++i) s[i] = double(v[i]) + d;
		return s;
	}
	inline DSUM operator-(double d) const
	{
		DSUM s;
		for (int i = 0; i < N; ++i) s[i] = double(v[i]) - d;
		return s;
	}
	inline DSUM operator*(double d) const
	{
		DSUM s;
		for (int i = 0; i < N; ++i) s[i] = double(v[i]) * d;
		return s;
	}
	inline DSUM operator/(double d) const
	{
		DSUM s;
		for (int i = 0; i < N; ++i) s[i] = double(v[i]) / d;
		return s;
	}
	inline ptvn& operator+=(const ptvn& t)
	{
		for (int i = 0; i < N; ++i) type_copy(v[i] + t[i], v[i]);
		return *this;
	}
	inline ptvn& operator-=(const ptvn& t)
	{
		for (int i = 0; i < N; ++i) type_copy(v[i] - t[i], v[i]);
		return *this;
	}
	inline ptvn& operator*=(const ptvn& t)
	{
		for (int i = 0; i < N; ++i) type_copy(v[i] * t[i], v[i]);
		return *this;
	}
	inline ptvn& operator/=(const ptvn& t)
	{
		for (int i = 0; i < N; ++i) type_copy(v[i] / t[i], v[i]);
		return *this;
	}
	inline ptvn& operator+=(double d)
	{
		for (int i = 0; i < N; ++i) type_copy(v[i] + d, v[i]);
		return *this;
	}
	inline ptvn& operator-=(double d)
	{
		for (int i = 0; i < N; ++i) type_copy(v[i] - d, v[i]);
		return *this;
	}
	inline ptvn& operator*=(double d)
	{
		for (int i = 0; i < N; ++i) type_copy(v[i] * d, v[i]);
		return *this;
	}
	inline ptvn& operator/=(double d)
	{
		for (int i = 0; i < N; ++i) type_copy(v[i] / d, v[i]);
		return *this;
	}
	inline ptvn operator&(const ptvn& t) const
	{
		ptvn pt;
		for (int i = 0; i < N; ++i) pt[i] = v[i] & t[i];
		return pt;
	}
	inline ptvn operator|(const ptvn& t) const
	{
		ptvn pt;
		for (int i = 0; i < N; ++i) pt[i] = v[i] | t[i];
		return pt;
	}
	inline ptvn operator~() const
	{
		ptvn pt;
		for (int i = 0; i < N; ++i) pt[i] = ~v[i];
		return pt;
	}
	inline ptvn& operator&=(const ptvn& t)
	{
		for (int i = 0; i < N; ++i) v[i] &= t[i];
		return *this;
	}
	inline ptvn& operator|=(const ptvn& t)
	{
		for (int i = 0; i < N; ++i) v[i] |= t[i];
		return *this;
	}
	inline ptvn& operator=(const DSUM& s)
	{
		for (int i = 0; i < N; ++i) type_copy(s[i], v[i]);
		return *this;
	}
	inline ptvn& operator=(const ptvn& t)
	{
		for (int i = 0; i < N; ++i) v[i] = t[i];
		return *this;
	}
	inline bool operator==(const ptvn& t) const
	{
		for (int i = 0; i < N; ++i) if (v[i] != t[i]) return false;
		return true;
	}
	inline bool operator<(const ptvn& t) const
	{
		for (int i = 0; i < N; ++i) if (v[i] >= t[i]) return false;
		return true;
	}
	inline bool operator<=(const ptvn& t) const
	{
		for (int i = 0; i < N; ++i) if (v[i] > t[i]) return false;
		return true;
	}
	inline bool operator>(const ptvn& t) const
	{
		for (int i = 0; i < N; ++i) if (v[i] <= t[i]) return false;
		return true;
	}
	inline bool operator>=(const ptvn& t) const
	{
		for (int i = 0; i < N; ++i) if (v[i] < t[i]) return false;
		return true;
	}
	inline bool operator==(const DSUM& t) const
	{
		for (int i = 0; i < N; ++i) if (v[i] != t[i]) return false;
		return true;
	}
	inline bool operator<(const DSUM& t) const
	{
		for (int i = 0; i < N; ++i) if (v[i] >= t[i]) return false;
		return true;
	}
	inline bool operator<=(const DSUM& t) const
	{
		for (int i = 0; i < N; ++i) if (v[i] > t[i]) return false;
		return true;
	}
	inline bool operator>(const DSUM& t) const
	{
		for (int i = 0; i < N; ++i) if (v[i] <= t[i]) return false;
		return true;
	}
	inline bool operator>=(const DSUM& t) const
	{
		for (int i = 0; i < N; ++i) if (v[i] < t[i]) return false;
		return true;
	}
	inline double distance(const ptvn& t) const
	{
		DSUM s;
		for (int i = 0; i < N; ++i) s[i] = double(v[i]) - t[i];
		return s.mod();
	}
	inline double distance2(const ptvn& t) const
	{
		DSUM s;
		for (int i = 0; i < N; ++i) s[i] = double(v[i]) - t[i];
		return s.mod_square();
	}
	inline double average() const
	{
		double s = 0.0;
		for (int i = 0; i < N; ++i) s += v[i];
		return s / N;

	}
	inline T vmin() const
	{
		for (int i = 0; i < N; ++i)
		{
			int j = 0;
			for (j = 0; j < N; ++j) if (v[i] > v[j]) break;
			if (j == N) return v[i];
		}
	}
	inline T vmax() const
	{
		for (int i = 0; i < N; ++i)
		{
			int j = 0;
			for (j = 0; j < N; ++j) if (v[i] < v[j]) break;
			if (j == N) return v[i];
		}
	}
	inline DSUM gamma(double c, double g) const
	{
		DSUM s;
		for (int i = 0; i < N; ++i) s[i] = c * pow(double(v[i]), g);
		return s;
	}
	inline DSUM clog(double c) const
	{
		DSUM s;
		for (int i = 0; i < N; ++i) s[i] = c * log(1.0 + v[i]);
		return s;
	}
	T v[N];
};


template<typename TV1, typename TV2> void tv_convert(const TV1& v1, TV2& v2)
{
	int c1 = TV1::channels(), c2 = TV2::channels();
	if (c1 == c2)
	{
		for (int c = 0; c < c1; ++c) type_copy(v1[c], v2[c]);
	}
	else if (c1 < c2)
	{
		if (c1 == 1)
		{
			for (int c = 0; c < c2; ++c) type_copy(v1[0], v2[c]);
		}
		else
		{
			for (int c = 0; c < c1; ++c) type_copy(v1[c], v2[c]);
		}
	}
	else
	{
		if (c2 == 1)
		{
			for (int c = 0; c < c1; ++c) type_copy(v1.average(), v2[0]);
		}
		else
		{
			for (int c = 0; c < c2; ++c) type_copy(v1[c], v2[c]);
		}
	}
}

template<typename TV1,typename TV2> void ptv_convert(const TV1* pin,TV2* pout,int len)
{
	int c1 = TV1::channels(), c2 = TV2::channels();
	if (c1 == c2)
	{
		for (int c = 0; c < c1; ++c)
		{
			for (int i = 0; i < len; ++i) type_copy(pin[i][c], pout[i][c]);
		}
	}
	else if (c1 < c2)
	{
		if (c1 == 1)
		{
			for (int c = 0; c < c2; ++c)
			{
				for (int i = 0; i < len; ++i) type_copy(pin[i][0], pout[i][c]);
			}
		}
		else
		{
			for (int c = 0; c < c1; ++c)
			{
				for (int i = 0; i < len; ++i) type_copy(pin[i][c], pout[i][c]);
			}
		}
	}
	else
	{
		if (c2 == 1)
		{
			for (int c = 0; c < c2; ++c)
			{
				for (int i = 0; i < len; ++i) type_copy(pin[i].average(), pout[i][0]);
			}
		}
		else
		{
			for (int c = 0; c < c2; ++c)
			{
				for (int i = 0; i < len; ++i) type_copy(pin[i][c], pout[i][c]);
			}
		}
	}
}

template<typename TV1, typename TV2> void ptv_channel_convert(const TV1* pin, int channel, TV2* pout, int len)
{
	int c2 = TV2::channels();
	for (int c = 0; c < c2; ++c)
	{
		for (int i = 0; i < len; ++i) type_copy(pin[i][channel], pout[i][c]);
	}
}

template<typename T> T img_mid_value(T* p, int len)
{
	for (int n = len; n > 1; --n)
	{
		for (int i = 1; i < n; ++i)
		{
			if (p[i] < p[i - 1])
			{
				T t = p[i];
				p[i] = p[i - 1];
				p[i - 1] = t;
			}
		}
	}
	return p[len / 2];
}

template<typename T> void gauss_channel(PTArea<T>& rect, int m)
{
	gauss_mask mask(m);
	int w = rect.width() - m + 1;
	int h = rect.height() - m + 1;
	double* pbuff = new double[w * m];
	double** pp = new double* [m];
	for (int i = 0; i < m; ++i) pp[i] = pbuff + (w * i);

	for (int y = 0; y < m - 1; ++y)
	{
		const T* p = rect.row(y);
		for (int x = 0; x < w; ++x)
		{
			double s=0.0;
			for (int j = 0; j < m; ++j) s += p[x + j] * mask[j];
			pp[y][x] = s;
		}
	}
	for (int y = 0; y < h; ++y)
	{
		const T* p = rect.row(y + m - 1);
		T* pt = rect.row(y + m / 2);
		for (int x = 0; x < w; ++x)
		{
			double s = 0.0;
			for (int j = 0; j < m; ++j) s += p[x + j] * mask[j];
			pp[m - 1][x] = s;
		}
		for (int x = 0; x < w; ++x)
		{
			double s = 0.0;
			for (int i = 0; i < m; ++i) s += pp[i][x] * mask[i];
			pt[x + m / 2] = T(s);
		}
		double* p0 = pp[0];
		for (int i = 0; i < m - 1; ++i) pp[i] = pp[i + 1];
		pp[m - 1] = p0;
	}
	delete[] pbuff;
	delete[] pp;
}

template<class T>
void g_d_channel(int wout, T* p1, T* p2, T* p3, PFIDual* pout)
{
	for (int x = 0; x < wout; ++x)
	{
		double gx=0.0, gy=0.0;
		gx -= p1[x];
		gy -= p1[x];
		gy -= p1[x + 1] * 2;
		gx += p1[x + 2];
		gy -= p1[x + 2];
		gx -= p2[x] * 2;
		gx += p2[x + 2] * 2;
		gx -= p3[x];
		gy += p3[x];
		gy += p3[x + 1] * 2;
		gx += p3[x + 2];
		gy += p3[x + 2];

		pout[x].v2 = 2;//记录梯度方向
		pout[x].v1 = sqrt(gx*gx+gy*gy);
		if (gx != 0.0) pout[x].v2 = pdsum_dirt[int(atan(gy / gx) * 8.0 / pi) + 4];
	}
}

template<class T>
void g_m_channel(const T line1[3], const T line2[3], const T line3[3], PDSum3& mt)
{
	double gx = 0.0, gy = 0.0;
	gx -= line1[0];
	gy -= line1[0];
	gy -= line1[1] * 2;
	gx += line1[2];
	gy -= line1[2];
	gx -= line2[0] * 2;
	gx += line2[2] * 2;
	gx -= line3[0];
	gy += line3[0];
	gy += line3[1] * 2;
	gx += line3[2];
	gy += line3[2];
	mt[0] = gx * gx;//x^2
	mt[1] = gy * gy;//y^2
	mt[2] = gx * gy;//x*y
}

template<class TV>
void gradient_direction(int wout,TV* p1,TV* p2,TV* p3,PFIDual* pout)
{
	for(int x=0;x<wout;++x)
	{
		TV::DSUM gx,gy;
		gx-=p1[x];
		gy-=p1[x];
		gy-=p1[x+1]*2;
		gx+=p1[x+2];
		gy-=p1[x+2];
		gx-=p2[x]*2;
		gx+=p2[x+2]*2;
		gx-=p3[x];
		gy+=p3[x];
		gy+=p3[x+1] * 2;
		gx+=p3[x+2];
		gy+=p3[x+2];
	
		pout[x]=gx.gradient_direction(gy);
	}
}

template<class TV>
void gradient_moment(const TV line1[3], const TV line2[3], const TV line3[3], PDSum3& mt)
{
	TV::DSUM gx, gy;
	gx -= line1[0];
	gy -= line1[0];
	gy -= line1[1] * 2;
	gx += line1[2];
	gy -= line1[2];
	gx -= line2[0] * 2;
	gx += line2[2] * 2;
	gx -= line3[0];
	gy += line3[0];
	gy += line3[1] * 2;
	gx += line3[2];
	gy += line3[2];
	mt[0] = (gx * gx).sum();//x^2
	mt[1] = (gy * gy).sum();//y^2
	mt[2] = (gx * gy).sum();//x*y
}

void canny_restrain(int wout, PFIDual* pv1, PFIDual* pv2, PFIDual* pv3, PFIDual* pout);

class PHistogramItem
{
public:
	PHistogramItem()
	{
		zero();
	}
	inline PHistogramItem & operator=(const PHistogramItem &phi)
	{
		x0=phi.x0;
		x1=phi.x1;
		for(int i=0;i<256;++i) p[i]=phi.p[i];
	}
	inline void zero()
	{
		x0=x1=0.0;
		for(int i=0;i<256;++i) p[i]=0.0;
	}
	inline double range() const
	{
		return x1-x0;
	}
	void getLines(PLinklist<PLine> &lines,double k=0.01) const;
	double getValue1(const PLinklist<PLine> &lines) const;
	double getValue2(const PLinklist<PLine> &lines) const;
	inline double getValue1() const
	{
		PLinklist<PLine> lines;
		getLines(lines);
		return getValue1(lines);
	}
	inline double getValue2() const
	{
		PLinklist<PLine> lines;
		getLines(lines);
		return getValue2(lines);
	}
	inline double y(double x)
	{
		int i = int((x - x0) / (x1 - x0) * 256);
		if (i < 0) return 0.0;
		else if (i >= 256) return p[255];
		else
		{
			PDPoint p1,p2(i,p[i]);
			if (i > 0) p1.set(i-1, p[i - 1]);
			PLine l(p1, p2);
			return l.intersection(PLine(x, 0, x, p[255])).y;
		}
	}
	double p[256];
	double x0,x1;
};

class PZoomItem
{
public:
	PZoomItem() :v1(0), t1(0.0), v2(0), t2(0.0)
	{
	}
	int v1, v2;
	double t1, t2;
};

class PZoomGroup
{
public:
	PZoomGroup(int cx, int cy)
	{
		rxt = new PZoomItem[cx];
		ryt = new PZoomItem[cy];
	}
	~PZoomGroup()
	{
		delete[] rxt;
		delete[] ryt;
	}
	void make(const PPoint& spoint, const PDPoint& ratio, int w, int h, const PRect& rect);
	inline int x1() const
	{
		return _rect.x1;
	}
	inline int x2() const
	{
		return _rect.x2;
	}
	inline int y1() const
	{
		return _rect.y1;
	}
	inline int y2() const
	{
		return _rect.y2;
	}
	const PZoomItem& rx(int x) const
	{
		return rxt[x];
	}
	const PZoomItem& ry(int y) const
	{
		return ryt[y];
	}
protected:
	PZoomItem* rxt;
	PZoomItem* ryt;
	PRect _rect;
};

#define PIMAGETYPE_U8 1//8位无符号
#define PIMAGETYPE_S8 2//8位有符号 1通道
#define PIMAGETYPE_U16 3//16位无符号 1通道
#define PIMAGETYPE_S16 4//16位有符号 1通道
#define PIMAGETYPE_U32 5//32位无符号 1通道
#define PIMAGETYPE_S32 6//32位有符号 1通道
#define PIMAGETYPE_F32 7//32位浮点 1通道
#define PIMAGETYPE_D64 8//64位浮点 1通道

#define PIMAGEFILE_BMP 1
#define PIMAGEFILE_PNG 2
#define PIMAGEFILE_TIFF 3
#define PIMAGEFILE_JPEG 4

typedef ptv1<byte> u8v1;
typedef ptv3<byte> u8v3;
typedef ptv4<byte> u8v4;

class PVImageLayer:public PList<PSegment>
{
public:
	PVImageLayer()
	{
		color = u8v3();
	}
	PVImageLayer(u8v3 c)
	{
		color = c;
	}
	u8v3 color;
};

class PVImage
{
public:
	//PVImage()
	//{
	//}
	~PVImage()
	{
		for(int i=0;i<_layers.count();++i) delete _layers[i];
	}
	inline int layers() const
	{
		return _layers.count();
	}
	inline PVImageLayer& layer(int l)
	{
		return *_layers[l];
	}
	inline const PVImageLayer& layer(int l) const
	{
		return *_layers[l];
	}
	void drawLine(const PLine& line, u8v3 color, double w);//w 线宽
	void drawCircle(const PCircle& circle, u8v3 color, double w = 0.0);
	void drawEllipse(const PEllipse2& ellipse, u8v3 color, double w = 0.0);
	void drawRectangle(const PRectangle& rect, u8v3 color, double w = 0.0);
	void drawPolygon(const PPointset& polygon, u8v3 color, double w = 0.0);
	void drawPolyline(const PPointset& polyline, u8v3 color, double w);
	void drawBezierline(const PSPline2& spline, u8v3 color, double w);
	void drawCloseBezier(const PCloseSPline2& spline, u8v3 color, double w = 0.0);
	bool drawFileFigure(const wchar_t* file);
protected:
	PList<PVImageLayer*> _layers;
};

class PImageInfoTemp
{
public:
	PImageInfoTemp()
	{
		width = 0;
		height = 0;
		channels = 0;
		type = 0;
		size = 0;
		data = 0;
	}
	inline void* row(int y)
	{
		return data + y * (width * channels * size);
	}
	inline const void* row(int y) const
	{
		return data + y * (width * channels * size);
	}
	byte* data;
	int width;
	int height;
	int channels;
	int type;
	int size;//每通道单元字节数
	int saveBmp(const TCHAR* fileName) const;
	int saveTiff(const TCHAR* fileName) const;
	int savePng(const TCHAR* fileName) const;
	int saveJpeg(const TCHAR* fileName) const;
};

class PImageInfo:public PImageInfoTemp
{
public:
	PImageInfo():PImageInfoTemp()
	{
	}
	~PImageInfo()
	{
		if(data) delete[] data;
	}
	inline void destory()
	{
		if (data) delete[] data;
		data = 0;
	}
	int readBmp(const TCHAR *fileName);
	int readTiff(const TCHAR *fileName);
	int readPng(const TCHAR *fileName);
	int readJpeg(const TCHAR *fileName);
};

class PImageBase
{
public:
	PImageBase()
	{
		_width = 0;
		_height = 0;
	}
	inline int width() const
	{
		return _width;
	}
	inline int height() const
	{
		return _height;
	}
protected:
	int _width;
	int _height;
};

template<typename TV> class PImage:public PImageBase
{
public:
	PImage()
	{
		_data=0;
		_width=0;
		_height=0;
		ref = new int;
		*ref = 0;
	}
	PImage(const PImage &img)
	{
		ref = img.ref;
		++(*ref);
		_data = img._data;
		_width = img.width();
		_height = img.height();
		//_data = 0;
		//_width = img.width();
		//_height = img.height();
		//int len = _width * _height;
		//_data = new TV[len];
		//for (int i = 0; i < len; ++i) _data[i] = img._data[i];
	}
	PImage(int w, int h)
	{
		ref = new int;
		*ref = 0;
		_width = w;
		_height = h;
		_data = new TV[w * h];
	}
	~PImage()
	{
		if (*ref == 0)
		{
			if (_data) delete[] _data;
			_data = 0;
			_width = 0;
			_height = 0;
			delete ref;
		}
		else --(*ref);
	}
	void reset(int w, int h)
	{
		if (*ref == 0)
		{
			int len = _width * _height;
			if (len != w * h)
			{
				if (_data) delete[] _data;
				_data = new TV[w * h];
			}
			_width = w;
			_height = h;
		}
		else
		{
			--(*ref);
			ref = new int;
			*ref = 0;
			_width = w;
			_height = h;
			_data = new TV[w * h];
		}
	}
	inline TV pixel(int x, int y) const
	{
		return _data[y * _width + x];
	}
	inline TV& pixel(int x, int y)
	{
		return _data[y * _width + x];
	}
	static inline int channels()
	{
		return TV::channels();
	}
	static inline size_t size()
	{
		return TV::size();
	}
	static inline const type_info& type()
	{
		return TV::type();
	}
	inline const TV* row(int y) const
	{
		return _data + y * _width;
	}
	inline TV* row(int y)
	{
		return _data + y * _width;
	}
	inline const TV* data() const
	{
		return _data;
	}
	inline TV* data()
	{
		return _data;
	}
	void copy(const PRect& rect, PImage& out) const
	{
		for (int y = rect.y1; y < rect.y2; ++y)
		{
			const TV* p = row(y);
			TV* pout = out.row(y);
			for (int x = rect.x1; x < rect.x2; ++x) pout[x] = p[x];
		}
	}
	void copy_channel(const PRect& rect, PImage& out,int channel) const
	{
		for (int y = rect.y1; y < rect.y2; ++y)
		{
			const TV* p = row(y);
			TV* pout = out.row(y);
			for (int x = rect.x1; x < rect.x2; ++x) pout[x][channel] = p[x][channel];
		}
	}
	PImage copy() const
	{
		PImage out(_width, _height);
		int len = _width * _height;
		for (int i = 0; i < len; ++i) out._data[i] = _data[i];
		return out;
	}
	template<typename TW> void convert(PImage<TW>& out) const
	{
		out.reset(_width, _height);
		ptv_convert(_data, out.data(), _width * _height);
	}
	template<typename TW> void channel_convert(int channel, PImage<TW>& out) const
	{
		out.reset(_width, _height);
		ptv_channel_convert(_data, channel, out.data(), _width * _height);
	}
	inline int area() const
	{
		return _height * _width;
	}
	void set(const TV* p, int len)
	{
		int len0 = _width * _height;
		len0 = len0 < len ? len0 : len;
		for (int i = 0; i < len0; ++i) _data[i] = p[i];
	}
	PImage putChannel(const PImage<ptv1<typename TV::ITEM>>& gray,int channel) const
	{
		PImage out(_width, _height);
		int w = _width < gray.width() ? _width : gray.width();
		int h = _height < gray.height() ? _height : gray.height();
		for (int y = 0; y < h; ++y)
		{
			const TV* p = row(y);
			TV* pt = out.row(y);
			const ptv1<typename TV::ITEM>* pi = gray.row(y);
			for (int x = 0; x < w; ++x)
			{
				pt[x] = p[x];
				pt[x][channel] = pi[x][0];
			}
		}
		return out;
	}
	void setChannel(const PImage<ptv1<typename TV::ITEM>>& gray, int channel)
	{
		int w = _width < gray.width() ? _width : gray.width();
		int h = _height < gray.height() ? _height : gray.height();
		for (int y = 0; y < h; ++y)
		{
			TV* p = row(y);
			const ptv1<typename TV::ITEM>* pi = gray.row(y);
			for (int x = 0; x < w; ++x)
			{
				p[x][channel] = pi[x][0];
			}
		}
	}
	//void setChannel(const typename TV::ITEM* p, int len, int channel)
	//{
	//	int len0 = _width * _height;
	//	len0 = len0 < len ? len0 : len;
	//	for (int i = 0; i < len0; ++i) _data[i][channel] = p[i];
	//}
	PImage<ptv1<typename TV::ITEM>> getChannel(int channel) const
	{
		PImage<ptv1<typename TV::ITEM>> out(_width, _height);
		int len = _width * _height;
		ptv1<typename TV::ITEM>* p = out.data();
		for (int i = 0; i < len; ++i) p[i][0] = _data[i][channel];
		return out;
	}
	PImage& operator=(const PImage& img)
	{
		if (*ref == 0)
		{
			delete[] _data;
			delete ref;
		}
		else --(*ref);
		ref = img.ref;
		++(*ref);
		_data = img._data;
		_width = img.width();
		_height = img.height();
		//if (this != &img)
		//{
		//	reset(img.width(), img.height());
		//	int len = _width * _height;
		//	for (int i = 0; i < len; ++i) _data[i] = img._data[i];
		//}
		return *this;
	}
	PImage & operator-=(const PImage &img)
	{
		int w = _width < img.width() ? _width : img.width();
		int h = _height < img.height() ? _height : img.height();
		for (int y = 0; y < h; ++y)
		{
			TV* p = row(y);
			TV* pi = img.row(y);
			for (int x = 0; x < w; ++x) p[x] -= pi[x];
		}
		return (*this);
	}
	PImage & operator+=(const PImage &img)
	{
		int w = _width < img.width() ? _width : img.width();
		int h = _height < img.height() ? _height : img.height();
		for (int y = 0; y < h; ++y)
		{
			TV* p = row(y);
			TV* pi = img.row(y);
			for (int x = 0; x < w; ++x) p[x] += pi[x];
		}
		return (*this);
	}
	PImage& operator*=(const PImage &img)
	{
		int w = _width < img.width() ? _width : img.width();
		int h = _height < img.height() ? _height : img.height();
		for (int y = 0; y < h; ++y)
		{
			TV* p = row(y);
			TV* pi = img.row(y);
			for (int x = 0; x < w; ++x) p[x] *= pi[x];
		}
		return (*this);
	}
	PImage& operator/=(const PImage &img)
	{
		int w = _width < img.width() ? _width : img.width();
		int h = _height < img.height() ? _height : img.height();
		for (int y = 0; y < h; ++y)
		{
			TV* p = row(y);
			TV* pi = img.row(y);
			for (int x = 0; x < w; ++x) p[x] /= pi[x];
		}
		return (*this);
	}
	PImage & operator&=(const PImage &img)
	{
		int w = _width < img.width() ? _width : img.width();
		int h = _height < img.height() ? _height : img.height();
		for (int y = 0; y < h; ++y)
		{
			TV* p = row(y);
			TV* pi = img.row(y);
			for (int x = 0; x < w; ++x) p[x] &= pi[x];
		}
		return (*this);
	}
	PImage & operator|=(const PImage &img)
	{
		int w = _width < img.width() ? _width : img.width();
		int h = _height < img.height() ? _height : img.height();
		for (int y = 0; y < h; ++y)
		{
			TV* p = row(y);
			TV* pi = img.row(y);
			for (int x = 0; x < w; ++x) p[x] |= pi[x];
		}
		return (*this);
	}
	PImage operator+(const PImage &img) const
	{
		PImage out(_width,_height);
		int w = _width < img.width() ? _width : img.width();
		int h = _height < img.height() ? _height : img.height();
		for (int y = 0; y < h; ++y)
		{
			TV* p = row(y);
			TV* pt = out.row(y);
			TV* pi = img.row(y);
			for (int x = 0; x < w; ++x) pt[x] = p[x] + pi[x];
		}
		return out;
	}
	PImage operator-(const PImage &img) const
	{
		PImage out(_width, _height);
		int w = _width < img.width() ? _width : img.width();
		int h = _height < img.height() ? _height : img.height();
		for (int y = 0; y < h; ++y)
		{
			TV* p = row(y);
			TV* pt = out.row(y);
			TV* pi = img.row(y);
			for (int x = 0; x < w; ++x) pt[x] = p[x] - pi[x];
		}
		return out;
	}
	PImage operator*(const PImage &img) const
	{
		PImage out(_width, _height);
		int w = _width < img.width() ? _width : img.width();
		int h = _height < img.height() ? _height : img.height();
		for (int y = 0; y < h; ++y)
		{
			TV* p = row(y);
			TV* pt = out.row(y);
			TV* pi = img.row(y);
			for (int x = 0; x < w; ++x) pt[x] = p[x] * pi[x];
		}
		return out;
	}
	PImage operator/(const PImage &img) const
	{
		PImage out(_width, _height);
		int w = _width < img.width() ? _width : img.width();
		int h = _height < img.height() ? _height : img.height();
		for (int y = 0; y < h; ++y)
		{
			TV* p = row(y);
			TV* pt = out.row(y);
			TV* pi = img.row(y);
			for (int x = 0; x < w; ++x) pt[x] = p[x] / pi[x];
		}
		return out;
	}
	PImage operator&(const PImage &img) const
	{
		PImage out(_width, _height);
		int w = _width < img.width() ? _width : img.width();
		int h = _height < img.height() ? _height : img.height();
		for (int y = 0; y < h; ++y)
		{
			TV* p = row(y);
			TV* pt = out.row(y);
			TV* pi = img.row(y);
			for (int x = 0; x < w; ++x) pt[x] = p[x] & pi[x];
		}
		return out;
	}
	PImage operator|(const PImage &img) const
	{
		PImage out(_width, _height);
		int w = _width < img.width() ? _width : img.width();
		int h = _height < img.height() ? _height : img.height();
		for (int y = 0; y < h; ++y)
		{
			TV* p = row(y);
			TV* pt = out.row(y);
			TV* pi = img.row(y);
			for (int x = 0; x < w; ++x) pt[x] = p[x] | pi[x];
		}
		return out;
	}
	PImage negation() const
	{
		PImage out(_width , _height);
		int len = _width * _height;
		for (int i = 0; i < len; ++i)
		{
			out._data[i] = ~_data[i];
		}
		return out;
	}
	PImage half() const
	{
		PImage out(_width/2, _height/2);
		for(int y=0;y<out.height();++y)
		{
			const TV* p0 = row(y * 2);
			const TV* p1 = row(y * 2 + 1);
			TV* pt=out.row(y);
			for(int x=0;x<out.width();++x)
			{
				TV::DSUM sum;
				sum += p0[x * 2];
				sum += p0[x * 2 + 1];
				sum += p1[x * 2];
				sum += p1[x * 2 + 1];
				pt[x]=sum*0.25;
			}
		}
		return out;
	}
	PImage zoom(int w,int h) const
	{
		PImage out(w,h);
		PDPoint ratio(double(w)/_width,double(h)/_height);
		zoom(PPoint(),ratio,out,PRect(0,0,w,h));
		return out;
	}
	void zoom(const PPoint& spoint,const PDPoint& ratio,PImage& out,const PRect& rect) const
	{
		double pix_area = ratio.area();
		PZoomGroup rt(out.width(),out.height());
		rt.make(spoint,ratio,width(),height(),rect);

		TV::DSUM* psum=new TV::DSUM[out.width()];
		for (int y = rt.y1(); y < rt.y2(); ++y)
		{
			const PZoomItem& ry = rt.ry(y);
			for (int x = rt.x1(); x < rt.x2(); ++x) psum[x].zero();
			const TV* pt = row(ry.v1);
			for (int x = rt.x1(); x < rt.x2(); ++x)
			{
				const PZoomItem& r = rt.rx(x);
				psum[x] += (pt[r.v1] * r.t1 + pt[r.v2] * r.t2) * ry.t1;
				for (int xt = r.v1 + 1; xt < r.v2; ++xt) psum[x] += pt[xt] * ry.t1;
			}
			for (int yt = ry.v1 + 1; yt < ry.v2; ++yt)
			{
				pt = row(yt);
				for (int x = rt.x1(); x < rt.x2(); ++x)
				{
					const PZoomItem& r = rt.rx(x);
					psum[x] += pt[r.v1] * r.t1 + pt[r.v2] * r.t2;
					for (int xt = r.v1 + 1; xt < r.v2; ++xt) psum[x] += pt[xt];
				}
			}
			pt = row(ry.v2);
			for (int x = rt.x1(); x < rt.x2(); ++x)
			{
				const PZoomItem& r = rt.rx(x);
				psum[x] += (pt[r.v1] * r.t1 + pt[r.v2] * r.t2) * ry.t2;
				for (int xt = r.v1 + 1; xt < r.v2; ++xt) psum[x] += pt[xt] * ry.t2;
			}
			TV* p = out.row(y);
			for (int x = rt.x1(); x < rt.x2(); ++x) p[x] = psum[x] * pix_area;//pix_area 缩放后的像素面积
		}
		delete[] psum;
	}
	void fastZoom(const PPoint& spoint,const PDPoint& ratio,PImage& out,const PRect& rect) const
	{
		int* ptemp = new int[out.width()];
		for (int x = rect.x1; x < rect.x2; ++x) ptemp[x] = int((x - spoint.x + 0.5) / ratio.x);

		int x1,x2;
		for (x1 = rect.x1; x1 < rect.x2; ++x1) if (ptemp[x1] >= 0) break;
		for (x2 = x1; x2 < rect.x2; ++x2) if (ptemp[x2] >= _width) break;

		for (int y = rect.y1; y < rect.y2; ++y)
		{
			TV* p = out.row(y);
			int yt = int((y - spoint.y + 0.5) / ratio.y);
			if (yt >= 0 && yt < _height)
			{
				for (int x = x1; x < x2; ++x)
				{
					p[x] = pixel(ptemp[x], yt);
				}
			}
		}
		delete[] ptemp;
	}
	PImage flip(int way) const//way 0 left-right,1 up-down,2 lu-rd,3 ru-ld
	{
		int w = width(), h = height();
		if(way==0)
		{
			PImage out(w, h);
			for(int y=0;y<h;++y)
			{
				const TV *p=row(y);
				TV *pt=out.row(y);
				for(int x=0;x<w;++x) pt[x]=p[w-1-x];
			}
			return out;
		}
		else if(way==1)
		{
			PImage out(w, h);
			for(int y=0;y<h;++y)
			{
				const TV *p=row(y);
				TV *pt=out.row(h-1-y);
				for(int x=0;x<w;++x) pt[x]=p[x];
			}
			return out;
		}
		else if(way==2)
		{
			PImage out(h, w);
			for(int y=0;y<w;++y)
			{
				TV *pt=out.row(y);
				for(int x=0;x<h;++x) pt[x]=pixel(w-1-y,h-1-x);
			}
			return out;
		}
		else //if(way==3)
		{
			PImage out(h, w);
			for(int y=0;y<w;++y)
			{
				TV *pt=out.row(y);
				for(int x=0;x<h;++x) pt[x]=pixel(y,x);
			}
			return out;
		}
	}
	void rotate(PImage& out,double sina,double cosa,const byte *ptab,const PRegion &reg) const
	{
		PDPoint o(_width*0.5,_height*0.5);
		PDPoint ot(out.width() * 0.5, out.height() * 0.5);
		PRegionPos pos(reg.begin());
		PRegionItem item;
		while(reg.getRange(pos,item))
		{
			TV *pt=out.row(item.y);
			PDPoint p(0.0,item.y+0.5-ot.y);
			for(int x=item.x0;x<item.x1;++x)
			{
				TV::DSUM s;
				p.x=x+0.5-ot.x;
				PDPoint q=p.rotate(-sina,cosa)+o;
				int x0=floor(q.x);//*floor
				int y0=floor(q.y);
		
				const TV *pn=row(y0-1);
				int x1=int((q.x-x0)*15.0);
				int y1=int((q.y-y0)*15.0);
				const byte *p9=ptab+y1*(15*9)+x1*9;
				s+=pn[x0-1]*p9[0];
				s+=pn[x0]*p9[1];
				s+=pn[x0+1]*p9[2];
				pn+=_width;
				s+=pn[x0-1]*p9[3];
				s+=pn[x0]*p9[4];
				s+=pn[x0+1]*p9[5];
				pn+=_width;
				s+=pn[x0-1]*p9[6];
				s+=pn[x0]*p9[7];
				s+=pn[x0+1]*p9[8];
				pt[x]=s/225.0;
			}
		}
	}
	PImage rotate(double angle,TV bk) const
	{
		PImage temp;
		PRegion reg;
		int wout,hout;
		expand(temp, 4, 4, bk);
		const byte *ptab=PImageRotatePara(_width,_height,angle,wout,hout,reg);
		PImage out(wout,hout);
		out.fill(bk);
		temp.rotate(out,sin(angle),cos(angle),ptab,reg);
		return out;
	}
	PImage affine(const matrix3x2 &mat) const//仿射变换
	{
	}
	PImage projection(const matrix3x3 &mat) const//投射变换
	{
	}
	PImage part(PRect rect) const
	{
		if(rect.x1>=_width || rect.y1>=_height) return;
		if(rect.x2>_width) rect.x2=_width;
		if(rect.y2>_height) rect.y2=_height;
		PImage out(rect.cx(),rect.cy());
		for(int y=0;y<rect.cy();++y)
		{
			const TV *p=row(y+rect.y1);
			TV *pout=row(y);
			for(int x=0;x<rect.cx();++x)
			{
				pout[x]=p[x+rect.x1];
			}
		}
		return out;
	}
	void fill(TV bk)
	{
		int len=area();
		for(int i=0;i<len;++i) _data[i]=bk;
	}
	void fill(const PRect &rect, TV bk)
	{
		PRect rt=rect;
		if(rt.x1>=_width || rt.y1>=_height) return;
		if(rt.x2>_width) rt.x2=_width;
		if(rt.y2>_height) rt.y2=_height;
		for(int y=rt.y1;y<rt.y2;++y)
		{
			TV* p=row(y)
			for(int x=rt.x1;x<rt.x2;++x) p[x]= bk;
		}
	}
	PImage gain(double mult,double add) const
	{
		PImage out(_width, _height);
		gain(PRect(0, 0, _width, _height), out, mult, add);
		return out;
	}
	void gain(const PRect &rect, PImage& out, double mult,double add) const
	{
		for(int y=rect.y1;y<rect.y2;++y)
		{
			const TV* p = row(y);
			TV *pt=out.row(y);
			for(int x=rect.x1;x<rect.x2;++x) pt[x]=p[x]*mult+add;
		}
	}
	PAreaItem areaCenter(const PRegion &region) const//返回0，表示未求出质心
	{
		PRegionPos pos(region.begin());
		PRegionItem item;
		double area=0.0;
		double wx=0.0,wy=0.0;
		while(region.getRange(pos,item))
		{
			double a=0.0;
			const TV *p=row(item.y);
			for(int x=item.x0;x<item.x1;++x)
			{
				wx+=(0.5+x)*p[x].average();
				a+=p[x];
			}
			wy+=a*(0.5+item.y);
			area+=a;
		}
		if(area>0.0) return PAreaItem(wx/area,wy/area,area);
		else return PAreaItem();
	}
	typename TV::DSUM meanValue() const//均值
	{
		TV::DSUM sum;
		for(int y=0;y< _height;++y)
		{
			TV::DSUM s=0.0;
			TV* p=row(y);
			for(int x=0;x< _widht;++x) s+=p[x];
			sum+=s/_widht;
		}
		return sum/h;
	}
	double maxValue() const//最大值
	{
		double v=-DBL_MAX;
		int len= area();
		for(int i=0;i<len;++i)
		{
			double a=_data[i].vmax();
			if(a>v) v=a;
		}
		return v;
	}
	double minValue() const//最小值
	{
		double v=DBL_MAX;
		int len= area();
		for(int i=0;i<len;++i)
		{
			double a=_data[i].vmin();
			if(a<v) v=a;
		}
		return v;
	}
	PDPair minmaxValue() const
	{
		PDPair p(DBL_MAX,-DBL_MAX);
		int len= area();
		for(int i=0;i<len;++i)
		{
			double a1=_data[i].vmin();
			double a2 = _data[i].vmax();
			if(a1<p.v1) p.v1=a1;
			if(a2>p.v2) p.v2=a2;
		}
		return p;
	}
	typename TV::DSUM stdDeviation(typename TV::DSUM vmean) const//标准差
	{
		TV::DSUM sum;
		for(int y=0;y<_height;++y)
		{
			TV::DSUM s=0.0;
			TV* p=row(y);
			for(int x=0;x<_widht;++x)
			{
				TV::DSUM m=vmean-p[x];
				s+=(vmean-p[x]()).square();
			}
			sum+=s/w;
		}
		return (sum/h).square_root();
	}
	PImage filter(int m,int n,const double *mask) const
	{
		PImage out(_width, _height);
		filter(PRect(0, 0, _width, _height), out, m, n, mask);
		return out;
	}
	void filter(const PRect& rect, PImage& out,int m,int n,const double *mask) const
	{
		PTArea<TV> temp;
		expand(rect,temp,m-1,n-1);
		int w = rect.cx();
		int h = rect.cy();
		for(int y=0;y<h;++y)
		{
			TV* p=temp.row(y);
			TV* pt = out.row(y+rect.y1);
			for(int x=0;x<w;++x)
			{
				TV::DSUM sum;
				const double *pm=mask;
				for(int i=0;i<n;++i)
				{
					TV *p1=temp.row(y+i)+x;
					for (int j = 0; j < m; ++j) sum += p1[j] * pm[j];
					pm+=n;
				}
				pt[x+rect.x1]=sum;
			}
		}
	}
	void gaussFilter(const PRect &rect,PImage& out,int m) const//多线程应用时,out不能为(*this)
	{
		gauss_mask mask(m);
		PTArea<TV> temp;
		expand(rect,temp,m-1,m-1);
		int w=rect.cx();
		int h=rect.cy();
		TV::DSUM* pbuff=new TV::DSUM[w*m];
		TV::DSUM** pp=new TV::DSUM*[m];
		for(int i=0;i<m;++i) pp[i]=pbuff+(w*i);

		for(int y=0;y<m-1;++y)
		{
			const TV* p=temp.row(y);
			for(int x=0;x<w;++x)
			{
				TV::DSUM s;
				for(int j=0;j<m;++j) s+=p[x+j]*mask[j];
				pp[y][x]=s;
			}
		}
		for(int y=0;y<h;++y)
		{
			const TV* p=temp.row(y+m-1);
			TV* pt=out.row(y+rect.y1)+rect.x1;
			for(int x=0;x<w;++x)
			{
				TV::DSUM s;
				for(int j=0;j<m;++j) s+=p[x+j]*mask[j];
				pp[m-1][x]=s;
			}
			for(int x=0;x<w;++x)
			{
				TV::DSUM s;
				for(int i=0;i<m;++i) s+=pp[i][x]*mask[i];
				pt[x]=s;
			}
			TV::DSUM* p0=pp[0];
			for(int i=0;i<m-1;++i) pp[i]=pp[i+1];
			pp[m-1]=p0;
		}
		delete[] pbuff;
		delete[] pp;
	}
	PImage gaussFilter(int m) const
	{
		PImage out(_width, _height);
		gaussFilter(PRect(0, 0, _width, _height), out, m);
		return out;
	}
	void meanSmooth(const PRect& rect,PImage& out,int m,int n) const
	{
		PTArea<TV> temp;
		expand(rect,temp,m-1,n-1);
		int w=rect.cx();
		int h=rect.cy();

		double mn=double(m*n);
		TV::DSUM* pbuff=new TV::DSUM[w*(n+1)];
		TV::DSUM** pp=new TV::DSUM*[n];
		for(int i=0;i<n;++i) pp[i]=pbuff+(w*i);
		TV::DSUM* psum=pbuff+(w*n);
		for(int y=0;y<n-1;++y)
		{
			TV *p=temp.row(y);
			TV::DSUM s;
			for(int j=0;j<m-1;++j) s+=p[j];
			for(int x=0;x<w;++x)
			{
				s+=p[x+m-1];
				pp[y][x]=s;
				s-=p[x];
			}
		}
		for(int x=0;x<w;++x)
		{
			TV::DSUM s;
			for(int i=0;i<n-1;++i) s+=pp[i][x];
			psum[x]=s;
		}
		for(int y=0;y<h;++y)
		{
			TV *p=temp.row(y+n-1);
			TV *pt=out.row(y+rect.y1) + rect.x1;
			TV::DSUM s;
			for(int j=0;j<m-1;++j) s+=p[j];
			for(int x=0;x<w;++x)
			{
				s+=p[x+m-1];
				pp[n-1][x]=s;
				s-=p[x];
			}
			for(int x=0;x<w;++x)
			{
				psum[x]+=pp[n-1][x];
				pt[x]=psum[x]/mn;
				psum[x]-=pp[0][x];
			}
			TV::DSUM* p0=pp[0];
			for(int i=0;i<n-1;++i) pp[i]=pp[i+1];
			pp[n-1]=p0;
		}
		delete[] pbuff;
		delete[] pp;
	}
	PImage meanSmooth(int m,int n) const
	{
		PImage out(_width, _height);
		meanSmooth(PRect(0, 0, _width, _height), out, m, n);
		return out;
	}
	void midSmooth(const PRect& rect,PImage& out,int m,int n) const
	{
		PTArea<TV> temp;
		expand(rect,temp,m-1,n-1);
		int w=rect.cx();
		int h=rect.cy();

		TV::DSUM* pbuff=new TV::DSUM[w*n];
		TV::DSUM** pp=new TV::DSUM*[n];
		TV::DSUM* pmid=new TV::DSUM[m+n];
		for(int i=0;i<n;++i) pp[i]=pbuff+(w*i);

		for(int y=0;y<n-1;++y)
		{
			TV* p=temp.row(y);
			for(int x=0;x<w;++x)
			{
				TV::DSUM s;
				for(int j=0;j<m;++j) pmid[j]=p[x+j];
				pp[y][x]=img_mid_value(pmid,m);
			}
		}
		for(int y=0;y<h;++y)
		{
			TV* p=temp.row(y+n-1);
			TV* pt=out.row(y+rect.y1)+rect.x1;
			for(int x=0;x<w;++x)
			{
				TV::DSUM s;
				for(int j=0;j<m;++j) pmid[j]=p[x+j];
				pp[n-1][x]=img_mid_value(pmid,m);
			}
			for(int x=0;x<w;++x)
			{
				TV::DSUM s;
				for(int i=0;i<n;++i) pmid[i]=pp[i][x];
				pt[x]=img_mid_value(pmid,n);
			}
			TV::DSUM* p0=pp[0];
			for(int i=0;i<n-1;++i) pp[i]=pp[i+1];
			pp[n-1]=p0;
		}
		delete[] pbuff;
		delete[] pp;
		delete[] pmid;
	}
	PImage midSmooth(int m,int n) const
	{
		PImage out(_width, _height);
		midSmooth(PRect(0, 0, _width, _height), out, m, n);
		return out;
	}
	PImage sobel() const
	{
		PImage out(_width, _height);
		sobel(PRect(0, 0, _width, _height), out);
		return out;
	}
	void sobel(const PRect& rect, PImage& out) const//rect计算区域,不检查越界
	{
		PTArea<TV> temp;
		expand(rect, temp, 2, 2);
		int w = rect.cx();
		int h = rect.cy();
		for (int y = 0; y < h; ++y)
		{
			const TV* p1 = temp.row(y);
			const TV* p2 = temp.row(y + 1);
			const TV* p3 = temp.row(y + 2);
			TV* pt = out.row(y + rect.y1) + rect.x1;
			for (int x = 0; x < w; ++x)
			{
				TV::DSUM gx, gy;
				gx -= p1[x];
				gy -= p1[x];
				gy -= p1[x + 1] * 2;
				gx += p1[x + 2];
				gy -= p1[x + 2];
				gx -= p2[x] * 2;
				gx += p2[x + 2] * 2;
				gx -= p3[x];
				gy += p3[x];
				gy += p3[x + 1] * 2;
				gx += p3[x + 2];
				gy += p3[x + 2];
				pt[x] = (gx.square() + gy.square()).square_root() * 0.25;
			}
		}
	}
	PImage laplace() const
	{
		double mask[9]={-1,-1,-1,-1,9,-1,-1,-1,-1};
		return filter(3,3,mask);
	}
	void canny(const PRect& rect,PMatRegion &out,double tlow,double thigh,int m) const
	{
		PTArea<TV> temp;
		expand(rect,temp,m+5,m+5);
		_gauss(temp,m);
		int w=rect.cx();
		int h=rect.cy();
		int w0=w+2; //h0=h+2;
		int w1=w+4; //h1=h+4;
		int w2=w+6; //h2=h+6;
		int x1=m/2;
		int y1=m/2;

		PFIDual *pu=new PFIDual[w0*3];
		PFIDual *pv=new PFIDual[w1*5];
		PFIDual *ppv[5],*ppu[3];
		for(int i=0;i<5;++i) ppv[i]=pv+(w1*i);
		for(int i=0;i<3;++i) ppu[i]=pu+(w0*i);

		for(int y=0;y<4;++y)
		{
			int yt=y+y1;
			gradient_direction(w1,temp.row(yt)+x1,temp.row(yt+1)+x1,temp.row(yt+2)+x1,ppv[y]);
		}
		canny_restrain(w0,ppv[0],ppv[1],ppv[2],ppu[0]);//非最大值抑制
		canny_restrain(w0,ppv[1],ppv[2],ppv[3],ppu[1]);

		for(int y=0;y<h;++y)
		{
			int yt=y+y1;
			byte *pout=out.row(y+rect.y1);
			gradient_direction(w1,temp.row(yt+4)+x1,temp.row(yt+5)+x1,temp.row(yt+6)+x1,ppv[4]);
			canny_restrain(w0,ppv[2],ppv[3],ppv[4],ppu[2]);//非最大值抑制

			PFIDual *pt=ppu[1];
			PFIDual *p1=ppu[0];
			PFIDual *p2=ppu[2];
			for(int x=0;x<w;++x)//滞后双阈值处理
			{
				int xt=x+rect.x1;
				if(pt[x+1].v1>thigh) pout[xt]=255;
				else if(pt[x+1].v1>=tlow)
				{
					 if(pt[x].v1>thigh || pt[x+2].v1>thigh
						|| p1[x].v1>thigh || p1[x+1].v1>thigh || p1[x+2].v1>thigh
						|| p2[x].v1>thigh || p2[x+1].v1>thigh || p2[x+2].v1>thigh) pout[xt]=255;
				}
				else pout[xt]=0;
			}
			pt=ppu[0];
			ppu[0]=ppu[1];
			ppu[1]=ppu[2];
			ppu[2]=pt;
			pt=ppv[0];
			for(int i=0;i<4;++i) ppv[i]=ppv[i+1];
			ppv[4]=pt;
		}
		delete[] pu;
		delete[] pv;
	}
	void canny(PMatRegion &out,double thigh,int m) const
	{
		double tlow=thigh*0.4;//确定低阈值
		if(m<0) m=(_width+_height)/300;//估计滤波模板
		m=(m/2)*2+1;//保证m是奇数
		if(m<3) m=3;
		out.reset(_width,_height);
		canny(PRect(0, 0, _width, _height), out, tlow, thigh, m);
		out.thinning();
	}
	PRegion canny(double thigh,int m) const
	{
		PMatRegion temp;
		canny(temp,thigh,m);
		return temp.threshold();
	}
	PRegion marrHildreth(double sigma,int m) const
	{
		PMatRegion temp(_width,_height);
		marrHildreth(temp, sigma, m);
		return temp.threshold();
	}
	void marrHildreth(PMatRegion &out,double sigma,int m) const
	{
		out.reset(_width,_height);
		marrHildreth(PRect(0, 0, _width, _height), out, sigma, m);
		out.thinning();
	}
	void marrHildreth(const PRect& rect, PMatRegion& out, double sigma, int m) const
	{
		PTArea<TV> temp;
		expand(rect, temp, m + 1, m + 1);
		int w1 = rect.width() + 2;
		//int h1 = rect.hieght() + 2;
		TV::DSUM* pi = new TV::DSUM[w1 * 3];
		TV::DSUM* pp[3] = { pi,pi + w1,pi + w1 * 2 };
		double* mask = new double[m * m];
		LoGMask(mask, sigma, m);

		for (int y = 0; y < 2; ++y)
		{
			for (int x = 0; x < w1; ++x)
			{
				TV::DSUM sum;
				for (int i = 0; i < m; ++i)
				{
					TV* p = temp.row(y + i) + x;
					for (int j = 0; j < m; ++j) sum += p[j] * mask[i * m + j];
				}
				pp[y][x] = sum;
			}
		}
		for (int y = 0; y < rect.height(); ++y)
		{
			byte* pt = out.row(y+rect.y1)+rect.x1;
			for (int x = 0; x < w1; ++x)
			{
				TV::DSUM sum;
				for (int i = 0; i < m; ++i)
				{
					TV* p = temp.row(y + i + 2) + x;
					for (int j = 0; j < m; ++j) sum += p[j] * mask[i * m + j];
				}
				pp[2][x] = sum;
			}
			TV::DSUM z;
			TV::DSUM* p0 = pp[0];
			TV::DSUM* p1 = pp[1];
			TV::DSUM* p2 = pp[2];
			for (int x = 0; x < rect.width(); ++x)
			{
				TV::DSUM nt1 = p0[x] + p0[x + 1] + p1[x] + p1[x + 1];
				TV::DSUM nt2 = p0[x + 1] + p0[x + 2] + p1[x + 1] + p1[x + 2];
				TV::DSUM nt3 = p1[x] + p1[x + 1] + p2[x] + p2[x + 1];
				TV::DSUM nt4 = p1[x + 1] + p1[x + 2] + p2[x + 1] + p2[x + 2];
				if (!(nt1 * nt2 >= z) || !(nt1 * nt3 >= z) || !(nt1 * nt4 >= z) || 
					!(nt2 * nt3 >= z) || !(nt2 * nt4 >= z) || !(nt3 * nt4 >= z)) pt[x] = 255;
				else pt[x] = 0;
			}
			pp[0] = pp[1];
			pp[1] = pp[2];
			pp[2] = p0;
		}
		delete[] mask;
		delete[] pi;
	}
	void imgVector(PTArea<PTPoint<char>>& out) const
	{
		out.reset(_width,_height);
		PTArea<TV> temp;
		expand(PRect(0,0, _width, _height),temp,2,2);
		for(int y=0;y<_height;++y)
		{
			TV *p1=temp.row(y);
			TV *p2=temp.row(y+1);
			TV *p3=temp.row(y+2);
			PTPoint<char> *pt=out.row(y);
			for(int x=0;x<_width;++x)
			{
				TV::DSUM vx,vy;
				vx+=p1[x];
				vx-=p1[x+2];
				vy+=p1[x];
				vy+=p1[x+1]*2;
				vy+=p1[x+2];
				vx+=p2[x]*2;
				vx-=p2[x+2]*2;
				vx+=p3[x];
				vx-=p3[x+2];
				vy-=p3[x];
				vy-=p3[x+1]*2;
				vy-=p3[x+2];
				PDPoint pv=vx.vector(vy);
				double m=pv.mod();
				if(m>0.0) pv/=m;				
				pt[x].set(char(pv.x*127.0),char(pv.y*127.0));
			}
		}
	}
	PDPoint pointVector(int x,int y) const
	{
		TV::DSUM vx,vy;
		const TV *p1=0;
		const TV *p2=0;
		const TV *p3=0;
		if((x==0 && y==0) || (x==_width-1 && y==_height-1)) return PDPoint();
		else if(x==0)
		{
			p1=row(y-1);
			p2=row(y);
			p3=row(y+1);
			vy+=p1[1];
			vy+=p1[0]*2;
			vy+=p1[1];
			vy-=p3[1];
			vy-=p3[0]*2;
			vy-=p3[1];
		}
		else if(x==_width-1)
		{
			p1=row(y-1)+x-1;
			p2=row(y)+x-1;
			p3=row(y+1)+x-1;
			vy+=p1[0];
			vy+=p1[1]*2;
			vy+=p1[0];
			vy-=p3[0];
			vy-=p3[1]*2;
			vy-=p3[0];
		}
		else if(y==0)
		{
			p1=row(0)+x-1;
			p2=row(1)+x-1;
			vx+=p2[0]*2;
			vx-=p2[2]*2;
			vx+=p1[0]*2;
			vx-=p1[2]*2;
		}
		else if(y==_height-1)
		{
			p1=row(y-1)+x-1;
			p2=p1+_width;
			vx+=p1[0]*2;
			vx-=p1[2]*2;
			vx+=p2[0]*2;
			vx-=p2[2]*2;
		}
		else
		{
			p1=row(y-1)+x-1;
			p2=p1+ _width;
			p3=p2+ _width;
			vx+=p1[0];
			vx-=p1[2];
			vy+=p1[0];
			vy+=p1[1]*2;
			vy+=p1[2];
			vx+=p2[0]*2;
			vx-=p2[2]*2;
			vx+=p3[0];
			vx-=p3[2];
			vy-=p3[0];
			vy-=p3[1]*2;
			vy-=p3[2];
		}
		PDPoint pv=vx.vector(vy);
		double m=pv.mod();
		if(m>0.0) pv/=m;
		return pv;
	}
	void harris(PPointset &pset,double thr,int m) const//角点检测,m 模板大小,thr 阈值
	{
		harris(PRect(0, 0, _width, _height), pset, thr, m);
	}
	void harris(const PRect& rect, PPointset& pset, double thr, int m) const//角点检测,m 模板大小,thr 阈值
	{
		pset.clear();
		PTArea<TV> temp;
		expand(rect, temp, m + m + 2, m + m + 2);
		_gauss(temp, m);
		int w = rect.cx();
		int h = rect.cy();
		int w1 = w + 2; 
		int w2 = w + m + 1;

		double* pr = new double[w1 * 3];
		double* ppr[3] = { pr,pr + w1,pr + w1 * 2 };
		PDSum3* pa = new PDSum3[w2 * (m + 1)];
		PDSum3** pp = new PDSum3*[m];
		for (int i = 0; i < m; ++i) pp[i] = pa + w2 * i;
		PDSum3* pb = pa + w2 * m;

		for (int y = 0; y < m; ++y)
		{
			TV* p1 = temp.row(y + m / 2) + m / 2;
			TV* p2 = temp.row(y + 1 + m / 2) + m / 2;
			TV* p3 = temp.row(y + 2 + m / 2) + m / 2;
			for (int x = 0; x < w2; ++x) gradient_moment(p1 + x, p2 + x, p3 + x, pp[y][x]);
		}
		for (int x = 0; x < w2; ++x)
		{
			for (int y = 0; y < m; ++y) pb[x] += pp[y][x];
		}
		for (int y = 0; y < 2; ++y)
		{
			PDSum3* p0 = pp[0];
			for (int x = 0; x < w2; ++x) pb[x] -= p0[x];
			for (int i = 0; i < m - 1; ++i) pp[i] = pp[i + 1];
			pp[m - 1] = p0;
			TV* p1 = temp.row(y + m + m / 2) + m / 2;
			TV* p2 = temp.row(y + m + 1 + m / 2) + m / 2;
			TV* p3 = temp.row(y + m + 2 + m / 2) + m / 2;
			for (int x = 0; x < w2; ++x) gradient_moment(p1 + x, p2 + x, p3 + x, p0[x]);
			for (int x = 0; x < w2; ++x) pb[x] += p0[x];

			double* prt = ppr[y];
			PDSum3 xy;
			for (int x = 0; x < m; ++x) xy += pb[x];
			prt[0] = (xy[0] * xy[1] - xy[2] * xy[2]) - 0.05 * (xy[0] + xy[1]) * (xy[0] + xy[1]);
			for (int x = 1; x < w1; ++x)
			{
				xy -= pb[x];
				xy += pb[x + m];
				prt[x] = (xy[0] * xy[1] - xy[2] * xy[2]) - 0.05 * (xy[0] + xy[1]) * (xy[0] + xy[1]);
			}
		}
		double R = thr * thr * thr * thr * m * m;
		for (int y = 0; y < h; ++y)
		{
			PDSum3* p0 = pp[0];
			for (int x = 0; x < w2; ++x) pb[x] -= p0[x];
			for (int i = 0; i < m - 1; ++i) pp[i] = pp[i + 1];
			pp[m - 1] = p0;

			TV* p1 = temp.row(y + m + m / 2) + m / 2;
			TV* p2 = temp.row(y + m + 1 + m / 2) + m / 2;
			TV* p3 = temp.row(y + m + 2 + m / 2) + m / 2;
			for (int x = 0; x < w2; ++x) gradient_moment(p1 + x, p2 + x, p3 + x, p0[x]);
			for (int x = 0; x < w2; ++x) pb[x] += p0[x];

			double* prt = ppr[2];
			PDSum3 xy;
			for (int x = 0; x < m; ++x) xy += pb[x];
			prt[0] = (xy[0] * xy[1] - xy[2] * xy[2]) - 0.05 * (xy[0] + xy[1]) * (xy[0] + xy[1]);
			for (int x = 1; x < w1; ++x)
			{
				xy -= pb[x];
				xy += pb[x + m];
				prt[x] = (xy[0] * xy[1] - xy[2] * xy[2]) - 0.05 * (xy[0] + xy[1]) * (xy[0] + xy[1]);
			}

			prt = ppr[1];
			double* pr1 = ppr[0];
			double* pr2 = ppr[2];
			for (int x = 0; x < w; ++x)
			{
				if (prt[x+1] > R && prt[x+1] > prt[x] && prt[x+1] > prt[x+2] &&
					prt[x+1] > pr1[x] && prt[x+1] > pr1[x+1] && prt[x+1] > pr1[x+2] &&
					prt[x+1] > pr2[x] && prt[x+1] > pr2[x+1] && prt[x+1] > pr2[x+2])
				{
					pset.append(PDPoint(0.5 + rect.x1 + x, 0.5 + rect.y1 + y));
				}
			}
			ppr[0] = ppr[1];
			ppr[1] = ppr[2];
			ppr[2] = pr1;
		}
		delete[] pa;
		delete[] pr;
		delete[] pp;
	}
	PImage hsi() const
	{
		PImage out(_width, _height);
		int len = area();
		for (int n = 0; n < len; ++n)
		{
			byte h, s, i;
			_data[n].hsi(h, s, i);
			out._data[n].set(h, s, i);
		}
		return out;
	}
	PImage rgb() const
	{
		PImage out(_width, _height);
		int len = area();
		for (int n = 0; n < len; ++n)
		{
			byte r, g, b;
			_data[n].rgb(r, g, b);
			out._data[n].set(b, g, r);
		}
		return out;
	}
	void gamma(const PRect &rect,PImage& out,double c,double g) const
	{
		for(int y=rect.y1;y<rect.y2;++y)
		{
			const TV* p=row(y);
			TV* pt=out.row(y);
			for(int x=rect.x1;x<rect.x2;++x)
			{
				pt[x]=p[x].gamma(c,g);
			}
		}
	}
	PImage gamma(double c,double g) const
	{
		PImage out(_width, _height);
		gamma(PRect(0, 0, _width, _height), out, c, g);
		return out;
	}
	void img_log(const PRect& rect, PImage& out, double c) const
	{
		for (int y = rect.y1; y < rect.y2; ++y)
		{
			const TV* p = row(y);
			TV* pt = out.row(y);
			for (int x = rect.x1; x < rect.x2; ++x)
			{
				pt[x] = p[x].clog(c);
			}
		}
	}
	PImage img_log(double c) const
	{
		PImage out(_width, _height);
		img_log(PRect(0, 0, _width, _height), out, c);
		return out;
	}
	typename TV::DSUM mean(const PRegion &region) const
	{
		TV::DSUM sum;
		int a=0;
		PRegionPos pos(region.begin());
		PRegionItem item;
		while(region.getRange(pos,item))
		{
			const TV* p=row(item.y);
			for(int x=item.x0;x<item.x1;++x)
			{
				sum+=p[x];
			}
			a+=item.length();
		}
		if(a>0) return sum/a;
		return sum;
	}
	typename TV::DSUM stdDeviation(const PRegion &region,typename TV::DSUM v_mean) const//已知均值求方差
	{
		TV::DSUM sum;
		int a=0;
		PRegionPos pos(region.begin());
		PRegionItem item;
		while(region.getRange(pos,item))
		{
			const TV* p=row(item.y);
			for(int x=item.x0;x<item.x1;++x)
			{
				TV::DSUM m=v_mean-p[x];
				sum+=m*m;
			}
			a+=item.length();
		}
		if(a>0) return (sum/a).square_root();
		return sum;
	}
	typename TV::DSUM graySum(const PRegion &region) const
	{
		TV::DSUM sum;
		PRegionPos pos(region.begin());
		PRegionItem item;
		while(region.getRange(pos,item))
		{
			const TV *p=row(item.y);
			for(int x=item.x0;x<item.x1;++x)
			{
				sum+=p[x];
			}
		}
		return sum;
	}
	void pointset(PPointset &pset,TV tmin,TV tmax) const
	{
		for(int y=0;y<_height;++y)
		{
			const TV* p=row(y);
			for(int x=0;x<_width;++x)
			{
				if(p[x]>=tmin && p[x]<tmax)
				{
					pset.append(PDPoint(x+0.5,y+0.5));
				}
			}
		}
	}
	PRegion dynamicThreshold(int m,typename TV::DSUM tmin,typename TV::DSUM tmax) const//m 模板大小
	{
		PRegion out;
		PImage temp=meanSmooth(m,m);
		for(int y=0;y<_height;++y)
		{
			int x0=0,x1=0;
			const TV* p1=row(y);
			const TV* p2=temp.row(y);
			for(int x=0;x<_width;++x)
			{
				TV::DSUM t=p1[x]-p2[x];
				if(t>=tmin && t<tmax)
				{
					if(x0==x1)
					{
						x0=x;
						x1=x;
					}
					++x1;
				}
				else if(x0!=x1)
				{
					out.append(y,x0,x1);
					x0=0;
					x1=0;
				}
			}
			if(x0!=x1) out.append(y,x0,x1);
		}
		return out;
	}
	PRegion dynamicThreshold(const PRect& rect, int m, typename TV::DSUM tmin, typename TV::DSUM tmax) const//m 模板大小
	{
		PRegion out;
		PTArea<TV> temp;
		expand(rect, temp, m -1, m -1);
		_gauss(temp, m);
		int w = rect.width(), h = rect.height();
		for (int y = 0; y < h; ++y)
		{
			int x0 = 0, x1 = 0;
			const TV* p1 = row(y+rect.y1)+rect.x1;
			const TV* p2 = temp.row(y+m/2)+m/2;
			for (int x = 0; x < w; ++x)
			{
				TV::DSUM t = p1[x] - p2[x];
				if (t >= tmin && t < tmax)
				{
					if (x0 == x1)
					{
						x0 = x;
						x1 = x;
					}
					++x1;
				}
				else if (x0 != x1)
				{
					out.append(y, x0, x1);
					x0 = 0;
					x1 = 0;
				}
			}
			if (x0 != x1) out.append(y+rect.y1, x0+rect.x1, x1+rect.x1);
		}
		return out;
	}
	PRegion histogramThreshold(int section) const//直方图自动阈值分割 sec 0暗,1中间,2亮
	{
		int c = channels();
		PHistogramItem* phi = new PHistogramItem[c];
		histogram(phi);
		TV::DSUM tmin, tmax;
		for (int i = 0; i < c; ++i)
		{
			if (section == 0)
			{
				tmin[i] = -DBL_MAX;
				tmax[i] = phi[i].getValue1();
			}
			else if (section == 1)
			{
				tmin[i] = phi[i].getValue1();
				tmax[i] = phi[i].getValue2();
			}
			else if (section == 2)
			{
				tmin[i] = phi[i].getValue2();
				tmax[i] = DBL_MAX;
			}
		}
		delete[] phi;
		return threshold(tmin, tmax);
	}
	PRegion automaticThreshold(int section) const//边缘增强的直方图自动阈值分割  sec 0暗,1中间,2亮
	{
		int c = channels();
		PHistogramItem* phi = new PHistogramItem[c];
		histogram(phi);
		TV::DSUM tmin, tmax;
		for (int i = 0; i < c; ++i)
		{
			if (section == 0)
			{
				tmin[i] = -DBL_MAX;
				tmax[i] = phi[i].getValue1();
			}
			else if (section == 1)
			{
				tmin[i] = phi[i].getValue1();
				tmax[i] = phi[i].getValue2();
			}
			else if (section == 2)
			{
				tmin[i] = phi[i].getValue2();
				tmax[i] = DBL_MAX;
			}
		}
		delete[] phi;
		return threshold(tmin, tmax);
	}
	PRegion optimalThreshold(bool light) const//最优值分割
	{
		PRegion out;
		return out;
		//if(_type==PByteImage) img_optimal_threshold(_width,_height,light,_data_u8,out);
	}
	PRegion threshold(typename TV::DSUM tmin,typename TV::DSUM tmax) const//>=tmin,<tmax
	{
		PRegion out;
		for(int y=0;y<_height;++y)
		{
			int x0=0,x1=0;
			const TV* p=row(y);
			for(int x=0;x<_width;++x)
			{
				if(p[x]>=tmin && p[x]<tmax)
				{
					if(x0==x1)
					{
						x0=x;
						x1=x;
					}
					++x1;
				}
				else if(x0!=x1)
				{
					out.append(y,x0,x1);
					x0=0;
					x1=0;
				}
			}
			if(x0!=x1) out.append(y,x0,x1);
		}
		return out;
	}
	PRegion regionThreshold(const PRegion &region,typename TV::DSUM tmin,typename TV::DSUM tmax) const//>=tmin,<tmax
	{
		PRegion out;
		PRegionPos pos(region.begin());
		PRegionItem item;
		while(region.getRange(pos,item))
		{
			int x0=0,x1=0;
			const TV *p=row(item.y);
			for(int x=item.x0;x<item.x1;++x)
			{
				if(p[x]>=tmin && p[x]<tmax)
				{
					if(x0==x1)
					{
						x0=x;
						x1=x;
					}
					++x1;
				}
				else if(x0!=x1)
				{
					out.append(item.y,x0,x1);
					x0=0;
					x1=0;
				}
			}
			if(x0!=x1) out.append(item.y,x0,x1);
		}
		return out;
	}
	PEllipse2 ellipticAxis(const PRegion &region) const
	{
		PEllipse2 ellipse;
		PRegionPos pos(region.begin());
		PRegionItem item;
		PAreaItem aitem = areaCenter(region);
		if (aitem.area > 0.0)
		{
			double u20 = 0.0, u02 = 0.0, u11 = 0.0;
			while (region.getRange(pos, item))
			{
				const TV* p = row(item.y);
				double my1 = (item.y + 0.5) - aitem.y;
				double my2 = my1 * my1;
				u20 += my2 * item.length();
				for (int x = item.x0; x < item.x1; ++x)
				{
					u02 += ((x + 0.5) - aitem.x) * ((x + 0.5) - aitem.x);
					u11 += ((x + 0.5) - aitem.x) * my1;
				}
			}
			u20 /= aitem.area;
			u02 /= aitem.area;
			u11 /= aitem.area;
			ellipse.x = aitem.x;
			ellipse.y = aitem.y;
			double us = sqrt((u20 - u02) * (u20 - u02) + 4.0 * u11 * u11);
			ellipse.a = sqrt(2.0 * (u20 + u02 + us));
			ellipse.b = sqrt(2.0 * (u20 + u02 - us));
			ellipse.alpha = 0.5 * atan2(2.0 * u11, u02 - u20);
		}
		return ellipse;
	}
	void histogram(PHistogramItem phi[]) const
	{
		PDPair p = minmaxValue();
		double spa = 256.0 / (p.v2 - p.v1);
		int len = area();
		double a1 = 1.0 / len;
		for (int c = 0; c < channels(); ++c)
		{
			PHistogramItem& ph = phi[c];
			ph.zero();
			for (int i = 0; i < len; ++i)
			{
				int t = int(spa * (_data[i][c] - p.v1));
				ph.p[t] += 1.0;
			}
			ph.x0 = p.v1;
			ph.x1 = p.v2;
			for (int i = 1; i < 256; ++i) ph.p[i] += ph.p[i - 1];
			for (int i = 0; i < 256; ++i) ph.p[i] *= a1;
		}
	}
	void histogramEdge(PHistogramItem phi[]) const//边缘响应直方图
	{
		PImage img;
		gaussFilter(img,5);
		img.sobel();
		img.histogram(phi);
	}
	//void histogram()//直方图均衡
	//{
	//	double pa[256];
	//	byte pi[256];
	//	histogram(pa);
	//	for(int i=0;i<256;++i) pi[i]=byte(pa[i]*255.99);
	//	lut(pi);
	//}
	void drawLine(const PLine& line, TV color, double w)//w 线宽
	{
		PVImage vimg;
		u8v3 c;
		tv_convert(color, c);
		vimg.drawLine(line, c, w);
		mix(vimg);
	}
	void drawCircle(const PCircle& circle, TV color, double w = 0.0)
	{
		PVImage vimg;
		u8v3 c;
		tv_convert(color, c);
		vimg.drawCircle(circle, c, w);
		mix(vimg);
	}
	void drawEllipse(const PEllipse2& ellipse, TV color, double w = 0.0)
	{
		PVImage vimg;
		u8v3 c;
		tv_convert(color, c);
		vimg.drawEllipse(ellipse, c, w);
		mix(vimg);
	}
	void drawRectangle(const PRectangle& rect, TV color, double w = 0.0)
	{
		PVImage vimg;
		u8v3 c;
		tv_convert(color, c);
		vimg.drawRectangle(rect, c, w);
		mix(vimg);
	}
	void drawPolygon(const PPointset& polygon, TV color, double w = 0.0)
	{
		PVImage vimg;
		u8v3 c;
		tv_convert(color, c);
		vimg.drawPolygon(polygon, c, w);
		mix(vimg);
	}
	void drawPolyline(const PPointset& polyline, TV color, double w)
	{
		PVImage vimg;
		u8v3 c;
		tv_convert(color, c);
		vimg.drawPolyline(polyline, c, w);
		mix(vimg);
	}
	void drawBezierline(const PSPline2& spline, TV color, double w)
	{
		PVImage vimg;
		u8v3 c;
		tv_convert(color, c);
		vimg.drawBezierline(spline, c, w);
		mix(vimg);
	}
	void drawCloseBezier(const PCloseSPline2& spline, TV color, double w = 0.0)
	{
		PVImage vimg;
		u8v3 c;
		tv_convert(color, c);
		vimg.drawCloseBezier(spline, c, w);
		mix(vimg);
	}
	bool drawFileFigure(const wchar_t* file)
	{
		PVImage vimg;
		if (vimg.drawFileFigure(file))
		{
			mix(vimg);
			return true;
		}
		return false;
	}
	void _info(PImageInfoTemp& info) const
	{
		info.width = _width;
		info.height = _height;
		info.channels = channels();
		info.size = size();
		const type_info& tinfo = type();
		if (tinfo == typeid(byte)) info.type = PIMAGETYPE_U8;
		else if (tinfo == typeid(char)) info.type = PIMAGETYPE_S8;
		else if (tinfo == typeid(uint16_t)) info.type = PIMAGETYPE_U16;
		else if (tinfo == typeid(int16_t)) info.type = PIMAGETYPE_S16;
		else if (tinfo == typeid(uint32_t)) info.type = PIMAGETYPE_U32;
		else if (tinfo == typeid(int32_t)) info.type = PIMAGETYPE_S32;
		else if (tinfo == typeid(float)) info.type = PIMAGETYPE_F32;
		else if (tinfo == typeid(double)) info.type = PIMAGETYPE_D64;
		else info.type = 0;
		info.data = (byte*)_data;
	}
	int load(const PImageInfo& info)
	{
		int ret = 0;
		reset(info.width, info.height);
		int len= _width* _height;
		if (channels() < info.channels) ret = 2;
		if (size() < info.size) ret = 1;
		if (info.type == PIMAGETYPE_U8)
		{
			if(info.channels ==1) ptv_convert((ptv1<byte>*)info.data, _data, len);
			else if (info.channels == 2) ptv_convert((ptv2<byte>*)info.data, _data, len);
			else if (info.channels == 3) ptv_convert((ptv3<byte>*)info.data, _data, len);
			else if (info.channels == 4) ptv_convert((ptv4<byte>*)info.data, _data, len);
		}
		else if (info.type == PIMAGETYPE_S8)
		{
			if (info.channels == 1) ptv_convert((ptv1<char>*)info.data, _data, len);
			else if (info.channels == 2) ptv_convert((ptv2<char>*)info.data, _data, len);
			else if (info.channels == 3) ptv_convert((ptv3<char>*)info.data, _data, len);
			else if (info.channels == 4) ptv_convert((ptv4<char>*)info.data, _data, len);
		}
		else if (info.type == PIMAGETYPE_U16)
		{
			if (info.channels == 1) ptv_convert((ptv1<uint16_t>*)info.data, _data, len);
			else if (info.channels == 2) ptv_convert((ptv2<uint16_t>*)info.data, _data, len);
			else if (info.channels == 3) ptv_convert((ptv3<uint16_t>*)info.data, _data, len);
			else if (info.channels == 4) ptv_convert((ptv4<uint16_t>*)info.data, _data, len);
		}
		else if (info.type == PIMAGETYPE_S16)
		{
			if (info.channels == 1) ptv_convert((ptv1<int16_t>*)info.data, _data, len);
			else if (info.channels == 2) ptv_convert((ptv2<int16_t>*)info.data, _data, len);
			else if (info.channels == 3) ptv_convert((ptv3<int16_t>*)info.data, _data, len);
			else if (info.channels == 4) ptv_convert((ptv4<int16_t>*)info.data, _data, len);
		}
		else if (info.type == PIMAGETYPE_U32)
		{
			if (info.channels == 1) ptv_convert((ptv1<uint32_t>*)info.data, _data, len);
			else if (info.channels == 2) ptv_convert((ptv2<uint32_t>*)info.data, _data, len);
			else if (info.channels == 3) ptv_convert((ptv3<uint32_t>*)info.data, _data, len);
			else if (info.channels == 4) ptv_convert((ptv4<uint32_t>*)info.data, _data, len);
		}
		else if (info.type == PIMAGETYPE_S32)
		{
			if (info.channels == 1) ptv_convert((ptv1<int32_t>*)info.data, _data, len);
			else if (info.channels == 2) ptv_convert((ptv2<int32_t>*)info.data, _data, len);
			else if (info.channels == 3) ptv_convert((ptv3<int32_t>*)info.data, _data, len);
			else if (info.channels == 4) ptv_convert((ptv4<int32_t>*)info.data, _data, len);
		}
		else if (info.type == PIMAGETYPE_F32)
		{
			if (info.channels == 1) ptv_convert((ptv1<float>*)info.data, _data, len);
			else if (info.channels == 2) ptv_convert((ptv2<float>*)info.data, _data, len);
			else if (info.channels == 3) ptv_convert((ptv3<float>*)info.data, _data, len);
			else if (info.channels == 4) ptv_convert((ptv4<float>*)info.data, _data, len);
		}
		else if (info.type == PIMAGETYPE_D64)
		{
			if (info.channels == 1) ptv_convert((ptv1<double>*)info.data, _data, len);
			else if (info.channels == 2) ptv_convert((ptv2<double>*)info.data, _data, len);
			else if (info.channels == 3) ptv_convert((ptv3<double>*)info.data, _data, len);
			else if (info.channels == 4) ptv_convert((ptv4<double>*)info.data, _data, len);
		}
		return ret;
	}
	void put(PImageInfo& info)
	{
		PImageInfoTemp temp;
		_info(temp);
		info.destory();
		info.width = temp.width;
		info.height = temp.height;
		info.channels = temp.channels;
		info.size = temp.size;
		info.type = temp.type;
		int len = _width * _height;
		info.data = new byte[len * channels * size];
		TV* p = (TV*)info.data;
		for (int i = 0; i < len; ++i) p[i] = _data[i];
	}
	int readFile(const TCHAR* fileName, int fileTyep)
	{
		int ret = 0;
		PImageInfo file;
		if (fileTyep == PIMAGEFILE_BMP) ret = file.readBmp(fileName);
		else if (fileTyep == PIMAGEFILE_PNG) ret = file.readPng(fileName);
		else if (fileTyep == PIMAGEFILE_TIFF) ret = file.readTiff(fileName);
		else if (fileTyep == PIMAGEFILE_JPEG) ret = file.readJpeg(fileName);
		else ret = -3;
		if (ret == 0) ret=load(file);
		return ret;
	}
	int saveFile(const TCHAR* fileName, int fileTyep) const
	{
		int ret = 0;
		PImageInfoTemp file;
		_info(file);
		if (fileTyep == PIMAGEFILE_BMP) ret = file.saveBmp(fileName);
		else if (fileTyep == PIMAGEFILE_PNG) ret = file.savePng(fileName);
		else if (fileTyep == PIMAGEFILE_TIFF) ret = file.saveTiff(fileName);
		else if (fileTyep == PIMAGEFILE_JPEG) ret = file.saveJpeg(fileName);
		else ret = -3;//参数错误
		return ret;
	}
	void expand(PImage& out, int mx, int my, TV bk) const//扩展和填充背景,mx, my : 宽度和高度扩展的增量
	{
		int ml = mx / 2, mr = mx - ml;
		int mt = my / 2, mb = my - mt;
		out.reset(_width + mx, _height + my);
		for (int i = 0; i < mt; ++i)
		{
			TV* pt = out.row(i);
			for (int j = 0; j < out.width(); ++j) pt[j] = bk;
		}
		for (int i = 0; i < _height; ++i)
		{
			const TV* p = row(i);
			TV* pt = out.row(i+mt);
			for (int j = 0; j < ml; ++j) pt[j] = bk;
			for (int j = 0; j < _width; ++j) pt[j+ml] = p[j];
			for (int j = 0; j < mr; ++j) pt[_width+ml+j] = bk;
		}
		for (int i = 0; i < mb; ++i)
		{
			TV* pt = out.row(_height+mt+i);
			for (int j = 0; j < out.width(); ++j) pt[j] = bk;
		}
	}
protected:
	void expand(const PRect& rect,PTArea<TV>& out,int mx,int my) const//rect:源图像中的范围; mx,my:宽度和高度扩展的增量
	{
		int ml=mx/2,mr=mx-ml;
		int mt=my/2,mb=my-mt;
		out.reset(rect.cx()+mx,rect.cy()+my);

		PRect rect1;
		rect1.x1=rect.x1-ml;
		rect1.x2=rect.x2+mr;
		rect1.y1=rect.y1-mt;
		rect1.y2=rect.y2+mb;

		int* posx=new int[out.width()];
		int* posy=new int[out.height()];
		for(int j=0;j<out.width();++j)
		{
			int jt=j+rect1.x1;
			if(jt<0) jt=-(jt+1);
			else if(jt>=_width) jt=_width*2-jt-1;
			posx[j]=jt;
		}
		for(int i=0;i<out.height();++i)
		{
			int it=i+rect1.y1;
			if(it<0) it=-(it+1);
			else if(it>=_height) it=_height*2-it-1;
			posy[i]=it;
		}
		for(int i=0;i<out.height();++i)
		{
			const TV* p=row(posy[i]);
			TV* pt=out.row(i);
			for(int j=0;j<out.width();++j) pt[j]=p[posx[j]];
		}
		delete[] posx;
		delete[] posy;
	}
	void expand(const PRect& rect,PTArea<typename TV::ITEM>& out,int channel,int mx,int my) const//rect:源图像中的范围; mx,my:宽度和高度扩展的增量
	{
		int ml=mx/2,mr=mx-ml;
		int mt=my/2,mb=my-mt;
		out.reset(rect.cx()+mx,rect.cy()+my);

		PRect rect1;
		rect1.x1=rect.x1-ml;
		rect1.x2=rect.x2+mr;
		rect1.y1=rect.y1-mt;
		rect1.y2=rect.y2+mb;

		int* posx=new int[out.width()];
		int* posy=new int[out.height()];
		for(int j=0;j<out.width();++j)
		{
			int jt=j+rect1.x1;
			if(jt<0) jt=-(jt+1);
			else if(jt>=_width) jt=_width*2-jt-1;
			posx[j]=jt;
		}
		for(int i=0;i<out.height();++i)
		{
			int it=i+rect1.y1;
			if(it<0) it=-(it+1);
			else if(it>=_height) it=_height*2-it-1;
			posy[i]=it;
		}
		for(int i=0;i<out.height();++i)
		{
			const TV* p=row(posy[i]);
			TV::ITEM* pt=out.row(i);
			for(int j=0;j<out.width();++j)
			{
				int jt=posx[j];
				pt[j]=p[jt][channel];
			}
		}
		delete[] posx;
		delete[] posy;
	}
	void _gauss(PTArea<TV>& rect,int m) const//
	{
		gauss_mask mask(m);
		int w=rect.width()-m+1;
		int h=rect.height()-m+1;
		TV::DSUM* pbuff=new TV::DSUM[w*m];
		TV::DSUM** pp=new TV::DSUM*[m];
		for(int i=0;i<m;++i) pp[i]=pbuff+(w*i);

		for(int y=0;y<m-1;++y)
		{
			const TV* p=rect.row(y);
			for(int x=0;x<w;++x)
			{
				TV::DSUM s;
				for(int j=0;j<m;++j) s+=p[x+j]*mask[j];
				pp[y][x]=s;
			}
		}
		for(int y=0;y<h;++y)
		{
			const TV* p=rect.row(y+m-1);
			TV* pt=rect.row(y+m/2);
			for(int x=0;x<w;++x)
			{
				TV::DSUM s;
				for(int j=0;j<m;++j) s+=p[x+j]*mask[j];
				pp[m-1][x]=s;
			}
			for(int x=0;x<w;++x)
			{
				TV::DSUM s;
				for(int i=0;i<m;++i) s+=pp[i][x]*mask[i];
				pt[x+m/2]=s;
			}
			TV::DSUM* p0=pp[0];
			for(int i=0;i<m-1;++i) pp[i]=pp[i+1];
			pp[m-1]=p0;
		}
		delete[] pbuff;
		delete[] pp;
	}
	void mix(const PVImage& vimg)
	{
		byte* mask = new byte[_width];
		for (int l = 0; l < vimg.layers(); ++l)
		{
			const PVImageLayer& layer = vimg.layer(l);
			double ymin = DBL_MAX;
			double ymax = 0.0;
			for (int i = 0; i < layer.count(); ++i)
			{
				if (layer[i].ymin < ymin) ymin = layer[i].ymin;
				if (layer[i].ymax > ymax) ymax = layer[i].ymax;
			}
			int y0 = int(ymin);
			int y1 = int(ceil(ymax));
			if (y0 < 0) y0 = 0;
			if (y1 > _height) y1 = _height;
			if (y0 >= _height || y1 <= 0) return;
			int ynum = y1 - y0;
			PList<double>* pbuff = new PList<double>[ynum * 8];
			for (int i = 0; i < layer.count(); ++i)
			{
				double yt0 = ceil(layer[i].ymin * 8.0) / 8.0;
				double yt1 = layer[i].ymax;
				if (yt0 < 0.0) yt0 = 0.0;
				if (yt1 > _height) yt1 = _height;
				for (double y = yt0; y < yt1; y += 0.125)
				{
					double x = layer[i].x(y);
					int j = int((y - y0) * 8.0);
					pbuff[j].append(x);
				}
			}
			TV t_color;
			tv_convert(layer.color, t_color);
			for (int i = 0; i < ynum; ++i)
			{
				for (int x = 0; x < _width; ++x) mask[x] = 0;
				for (int j = 0; j < 8; ++j)
				{
					PList<double>& pt = pbuff[i * 8 + j];
					if (pt.count() >= 2)
					{
						pt.ascend();
						for (int k = 0; k < pt.count(); k += 2)
						{
							double d1 = pt[k], d2 = pt[k + 1];
							if (d2 > 0.0 && d1 < _width)
							{
								if (d1 < 0.0) d1 = 0.0;
								if (d2 >= _width) d2 = _width - 0.0625;
								int x1 = int(d1), x2 = int(d2);
								if (x1 == x2) mask[x1] += byte((d2 - d1) * 8.0);
								else
								{
									mask[x1] += byte((x1 + 1.0 - d1) * 8.0);
									mask[x2] += byte((d2 - x2) * 8.0);
									for (int x = x1 + 1; x < x2; ++x) mask[x] += 8;
								}
							}
						}
					}
				}
				TV* p = row(i + y0);
				for (int x = 0; x < _width; ++x)
				{
					p[x] = (t_color * mask[x] + p[x] * (64.0 - mask[x])) / 64.0;
				}
			}
			delete[] pbuff;
		}
		delete[] mask;
	}
	TV* _data;
	int* ref;
};
